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Chapter 1 


Introduction 


Inelastic behavior of materials is a subject that has fascinated engineers and researchers over 
many decades. Sustained interest in this topic is due to its direct relevance in numerous 
engineering applications. The real challenge is to accurately predict inelastic deformation 
when the material is subjected to thermal transients and complex mechanical loading. Gas 
turbine engines, power generation systems, and automotive engines are just a few examples 
of where materials experience complex thermo-mechanical loading. In these situations, 
the material is subjected to multiaxial states of stress, different loading rates and variable 
temperature conditions. The accurate representation of stress, strain and temperature fields 
within the structural components depends strongly on the mathematical representations or 
constitutive equations of the inelastic behavior of these materials at various temperature 
levels. Further, to be generally applicable, the constitutive equations must be expressed in 
tensor ial form. 

The total strain in the material can be separated into elastic and inelastic parts. Clas- 
sical decomposition of the inelastic component into time independent (plastic) and time 
dependent (creep) contributions facilitates modeling of material behavior. The plastic com- 
ponent is considered to be history dependent, i.e., the current plastic strain depends on the 
loading history. The creep component is treated as time dependent and is responsible for 
deformation when the external loads and the temperature are held constant. 

At elevated temperatures, however, the classical decomposition into plastic and creep strains 
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breaks down. In structural alloys creep and plastic strains occur simultaneously and inter- 
actively. This suggests that we need a unified constitutive model to represent deformation 
of metals at elevated temperatures. Many such models have been developed in the last few 
decades (e.g. Walker [1981], Bodner [1987], Helling and Miller [1987], Robinson and Duffy 
[1990], Freed and Walker [1995], Arnold et al [1996]). Most constitutive theories use internal 
variables and appropriate evolutionary equations for simultaneously treating all aspects of 
inelastic deformation including plasticity, creep and stress relaxation. These theories are 
based on thermodynamics of irreversible processes and will be elaborated on in Chapter 3. 


1.1 Materials for high temperature applications 

Safety is a primary concern in aerospace applications. Thus, all substructures must remain 
structurally sound during all phases of flight. For gas turbine engine components, materials 
having good strength properties over a wide range of temperatures are required. Metal- 
matrix composites (MMCs) and superalloys possess this quality and hence are suitable 
for this application. Owing to its complexity, the fabrication of MMCs is expensive and 
time consuming. Hence, superalloys like Inconel 718 (a Nickel-based alloy) and Haynes 188 
(a cobalt-based alloy) are more popular than MMCs in the aerospace industry for engine 
applications. Aged Inconel 718 is strengthened by precipitate hardening, in which fine 
particles of a second phase are dispersed throughout the grains, thereby slowing down the 
movement of dislocations. The physical mechanisms responsible for strengthening of this 
alloy are explained in Chapter 2. 


1.2 Objective of the present research 

Inconel 718 is a wrought nickel-based alloy that is precipitation hardenable to obtain high 
strength and is the material to be investigated in the present research. 

Recent experiments on Inconel 718 indicate that the yield stress and subsequent flow stress 
are significantly higher in compression than in tension (Gil et al [1999b]). This phenomenon, 
known as the strength differential (SD), is observed for a wide range of temperatures as 
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shown in Figure 1.1. Figure 1.1(a) shows higher yield and post-yield stresses in compression 
at 650°C. Yield loci at 23 °C and 650°C, shown in Figure 1.1(b), dearly show eccentricity 
in the compressive stress direction. The SD effect can be mathematically defined as 

SD = 2 ( — x ioo (1.1) 

+ OtJ 

where a c and at are yield strengths in uniaxial compression and uniaxial tension, respec- 
tively. At the initiation of yielding the SD has a high value; it then decreases, and finally 
reaches a constant value at higher inelastic strains. 

A yield function, based solely on the second invariant of deviatoric stress (J 2 ) is not able 
to correctly represent a material with a SD. Thus, the goal of the current work is to 
generalize a viscoplastic model, currently based solely on to account for the 
strength differential effect. 


1.3 Generalization of Plasticity Theory 


Having established the need to generalize the yield function to capture the SD effect, it 
is worthwhile to look at some yield functions and flow rules used by other researchers, 
that deviate from J 2 plasticity theory. This section starts with some basic elements of J 2 
plasticity followed by instances where deviations from the classical approach axe required 
to capture material behavior. 

1.3.1 Metal Plasticity 

(Hill [1950], Mendelson [1983], Lubliner [1990], Chen [1994], Khan and Huang [1995]) 

At the continuum level, the initial material behavior is taken to be linear elastic up to 
the yield point. The yield surface delineates the current elastic region in stress space. For 
an elastically isotropic metal, yielding depends on stress, temperature, and internal state 
(described by hardening or internal state variables). Under isothermal conditions, initial 
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(a) 



(b) 



Figure 1.1: Strength Differential in aged Inconel 718: (a) Uniaxial Response at 650° C (b) 
Yield Loci for a 30/xm/m offset strain definition 
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yielding depends only on the stress state, which is often described by three stress invariants: 


h = (?ii 

J2 — 

J3 — g S%j Sj k S^i 

where I\ is the first stress invariant, J 2 and J 3 are the second and third deviatoric stress 
invariants respectively, a ij and Sij denote the Cauchy and deviatoric stress components 
respectively. The yield function can be expressed as 


f(<Tij)=f(IuJ2,J3)=0. (1-2) 

Physically, I\ represents the hydrostatic stress and J 2 represents the distortional energy in 
the material. Though no definite physical quantity is attributed to J 3 , it can be treated as 
a weighting parameter that induces asymmetry in yield and flow behavior between tension 
and compression. 


Based on the work of Bridgman [1952], yielding in most metals (e.g. mild steel, copper, 
aluminum) was found to be relatively insensitive to hydrostatic pressure. This suggests 
that / 1 is unimportant in the definition of yield. Also, for many metals the tensile and 
compressive flow behavior are identical (e.g. aluminum, titanium, unaged maraging steel, 
unaged Inconel 718), making J 3 unimportant. Hence the yield function depends only on J 2 
and the yield surface is an infinite cylinder in the principal stress space with its axis along 
the hydrostatic stress line. The von Mises yield criterion takes the form of a cylinder with 
a circular cross-section, resulting in 

y/5, T 2 = (Ty (1.3) 

where a y is the uniaxial tensile yield strength. 


When applied to these materials, Druckers stability postulate (Drucker [1951] and [1959]) 
results in the inelastic strain rate vector being normal to the yield surface. Thus, the yield 
function can be treated as a plastic potential and the inelastic strain rate can be determined 
by use of an associated flow law, 


3 


df_ 

V 




(1.4) 
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where A is a scalar that is determined from the consistency condition discussed in Chapter 
3. Hence classical plasticity theory, when applied to metals, suggests that the yield function 
is based only on shear stress and that the flow law is associative. 

1.3.2 Deviations from classical plasticity 

Though many metals follow the von Mises (J 2 based) yield criterion and associative flow 
law, there are some that deviate from either or both. In this section, we present some 
instances where these deviations are observed. By doing so, we further emphasize the need 
for the proposed investigation. 

Starting from the middle 1960s, experimental results were published which show that yield 
and flow stresses in tension and compression are different for some high strength steels (e.g.. 
Leslie and Sober [1967]; Kalish and Cohen [1969]; Rauch and Leslie [1972]; Chait [1972]; 
Spitzig et al [1975]). The physical mechanisms that can cause a SD were reviewed by Hirth 
and Cohen [1970] and Drucker [1973]. For an elastically isotropic material exhibiting SD, 
a yield function that only depends on J 2 will not give an accurate representation of the 
material behavior. 

Experiments on iron single crystals (Spitzig [1979]) showed that contrary to Schmid’s law 
(Schmid [1924]), the normal stress on the slip planes affects the dislocation interactions. 
The results suggest that the predominant effect of hydrostatic pressure is to retard the 
generation of mobile dislocations, resulting in an increase in dislocation density which in 
turn causes an increase in volume. This suggests that application of hydrostatic pressure 
affects inelastic behavior of the material, indicating that flow is dependent on I\. 

The above examples are a few instances where a J 2 based yield description falls short in 
representing the material behavior. Chapter 3 reviews more instances where a generalization 
of the yield function is needed. 
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1.3.3 Modified Yield Functions 


Some researchers have extended J 2 plasticity theory to account for deviations from classical 
plasticity. Drucker [1949] compared Osgood’s [1947] experimental data on aluminum alloy 
tubes to classical Tresca and von Mises yield criterion. While neither the Tresca criterion, 
nor the Mises criterion, agreed well with the experimental data, a function of the form 

= 4 - 2-25 Ji - 1, (1.5) 

did fit the data well. This was one of the first functions to go beyond the J 2 representation 
of yield for metallic materials. 

For materials exhibiting SD, J 2 plasticity theory needs to be generalized. Yield functions 
for such materials must be represented by an odd power of 7 1? an odd power of J3, or a 
combination of the two. Dependence on I\ requires pressure sensitivity, while dependence 
on J3 (and not /1) requires pressure insensitivity. Spitzig et al [1975] superimposed tensile 
and compressive loads on hydrostatic pressure for different materials. The results clearly 
showed that the flow stress is pressure dependent. In many cases, the flow stress appears 
to be linearly related to hydrostatic pressure, suggesting the use of a Drucker-Prager yield 
criterion (Drucker and Prager [1952]); 

/ = ot I\ + yj~3 2 — k (1.6) 

where a and k are experimentally determined constants. In other cases, the relationship 
was not truly linear suggesting that either a J3 term or a different function of is required. 

Spitzig et al [1975] had good success using a Drucker-Prager type yield criterion of the form: 

\/3J2 = c — a I\ (1.7) 

where a and c are both strain dependent, but their ratio a/c(= a) is shown to be independent 
of strain. The SD predicted by equation 1.7 is simply 2a. 

Use of an associated flow law with a pressure dependent yield function results in a predicted 
permanent volume change approximately 15 times greater than that determined by pre-test 
and post-test density measurements and strain gage measurements (Spitzig et al [1975,1976]; 
Spitzig and Richmond [1980], Richmond and Spitzig [1984]). This suggests the use of a non- 
associated flow law with pressure dependent yield representations for these materials. 
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1.4 Method of Generalization 


The generalization of classical plasticity theory lies in not ignoring the dependence of flow 
on the stress invariants I\ and J 3 . This necessarily complicates the mathematics and raises 
the issue of what functional form of the stress invariants best represents inelastic flow. 

A general approach for determining the inelastic flow dependence on each of the three stress 
invariants would be to follow stress paths where only one of the three invariants varies. Two 
classical experiments that do this are: a hydrostatic pressure test, where only 7i varies; and a 
pure torsion test, where only J 2 varies. There are many other stress paths that are possible. 
Unfortunately, these stress paths generally require three-dimensional stress states, which 
are difficult to obtain in the laboratory. For experimental expediency, tests involving axial- 
torsional loading of thin-walled tubes are to be used in determining the significance of the 
first and third stress invariants, I\ and J 3 , respectively. 

Once the form of the threshold function is determined for IN718, characterization is done to 
quantify the material parameters in the model. Characterization experiments include shear 
tests and uniaxial (tensile and compressive) tests. Multiple combinations of stress invariants 
in the threshold function can result in an equally good correlation with the uniaxial test data. 
Experimental validation is required to help choose the proper combination of invariants. 
This is done by following load paths that extremize the differences between the predicted 
responses, based on the chosen threshold functions and by comparing the experiments with 
predictions. This approach can be applied to any material exhibiting a SD, e.g. martensitic 
steels, and intermetallics like titanium aluminide. 


1.5 Overview 


Various mechanisms of deformation in IN718 are discussed in Chapter 2. Also reviewed are 
the strengthening mechanisms in metallic alloys and the mechanisms responsible for SD in 
metals. By comparing and contrasting different mechanisms, an attempt is made to identify 
the ones that could be responsible for the observed SD effect in IN718. 

Theoretical background, explaining the fundamentals of rate-independent and rate-dependent 


N AS A/CR— 200 1-210715 


8 



plasticity with emphasis on metal plasticity, is presented in Chapter 3. Thermodynamic 
considerations leading to the use of internal state variables to account for various aspects 
of time dependent inelastic deformation are explained therein. Related work on various 
aspects of inelastic deformation, viz. yield and flow is also reviewed. 

The viscoplastic model to be used for this investigation is proposed in Chapter 4. Start- 
ing with the proposed threshold function, the flow and evolution laws are derived. The 
requirement of convexity of the yield function is discussed. The experimental program for 
determination of the material parameters in the proposed model is discussed in Chapter 5. 
A brief review of experimental work in metal plasticity is followed by a description of the 
experimental setup to be used in the present research. Test matrices for characterization 
and validation experiments are developed. 

Experimental results and the correlation between model and experiments are provided in 
Chapter 6. Six material parameter sets that fit the characterization data are presented. 
Choice of the right parameter set for Inconel 718 by correlating the predicted responses to 
biaxial validation tests and hydrostatic pressure tests is discussed in detail. This is followed 
by concluding remarks and suggested future work in Chapter 7. 
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Chapter 2 


Microstructure and deformation 
Mechanisms 


In this chapter, a detailed review of the microstructure and strengthening mechanisms in 
Inconel 718 is presented. Also documented are the mechanisms responsible for the SD 
effect in metallic materials in general. The effect of these mechanisms on the macroscopic 
behavior of metals is reviewed. Finally, by comparing and contrasting various mechanisms, 
an attempt is made to identify the ones responsible for the SD in Inconel 718. 

Numerous structural studies on Inconel 718 have been carried out in an attempt to correlate 
its properties to microstructure and heat treatment. This alloy has a large number of 
phases that form in a specific temperature range and have a characteristic morphology. 
The following is a list of different phases that are formed during various stages of heat 
treatment: 

1. MC phase: discrete particle phase formed on solidification and is stable up to 1200°C. 

2. Laves phase : round, island shaped particles that form on solidification in the high- 
niobium areas and are stable up to 1175°C. 

3. S phase: needle- like/plate-like that form on cooling during solidification in high nio- 
bium areas and are stable up to 1010°C in wrought alloys, and up to 1120°C in cast 
alloys, precipitating from 840° C to 980° C. 
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Table 2. 


: Material Composition of Inconel 718 


Element 

Content (wt.%) 

Ni 

53.58 

Cr 

17.52 

Mo 

2.87 

(Nb+Ta) 

5.19 

Ti 

0.95 

A1 

0.57 

Co 

0.39 

C 

0.034 

s 

0.002 

Mn 

0.12 

Si 

0.07 

B 

0.004 

Cu 

0.05 

P 

0.006 

Fe 

Bal. 


4. 7" phase: disc shaped precipitates form during cooling or heat treatment between 
730°C and 900°C. During cooling these precipitates tend to become smaller. 

5. 7* phase: spheroidal precipitates form when cooling or heat treatment is in the range, 
720°C-620°C. 

2.1 Microstructure of Inconel 718 

Wrought Inconel 718 alloy was obtained in the form of extruded 31.8 mm bars, all from 
the same heat treatment. The weight composition of the alloy is listed in Table 2.1. The 
machined samples were solutioned at 1038° C in argon for one hour and then air cooled. 
Then they were aged at 720°C in argon for eight hours, cooled at 55°C/hour to 620°C and 
held for eight hours, then air cooled to room temperature. 
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Figure 2.1: Optical Micrograph of aged Inconel 718 


Gil et al [1999b] performed metallography on this alloy by polishing and etching using 
Tucker’s reagent (45 m] HC1, 25 ml H 2 O, 15 ml HNO 3 , and 15 ml HF) to reveal the grain 
structure as seen in Figure 2.1. The grain structure consisted of equiaxed grains having 
an ASTM size of 4 (90 /im in diameter). Carbide particles were observed throughout the 
microstructure. Microhardness was measured to be Vicker 440 / Rockwell C of 45. 

Transmission electron microscopy showed a fine dispersion of 7 " precipitates with preferred 
orientation within a particular grain. This is shown in Figure 2.2 (samples prepared by C. 
Gil, observed by P. Howell, Earth and Mineral Sciences, PSU). The precipitate particles 
were observed to be platelets approximately 10-15 nm in length. Texture analysis (by M. 
Angelone, Material Characterization Lab, PSU) showed that there is no preferred grain ori- 
entation. Thus the preferred orientation of the precipitates within the grain is not expected 
to cause anisotropic behavior. 


2.2 Mechanisms of deformation in metals 

This section reviews some of the physical mechanisms that occur during inelastic deforma- 
tion in metallic materials. Strengthening mechanisms in aged Inconel 718 and the interac- 
tion between the applied stress and microstructure are discussed next. 
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Figure 2.2: Transmission Electron Microscopy of aged Inconel 718 showing 7 " precipitates 


At low temperature, slip is the most common mechanism of deformation. Slip occurs on 
the plane of highest atomic density in the direction closest to the planes of maximum shear 
stress. These planes and directions differ depending on the crystallographic structure of the 
metal and the direction of the applied load. 

During inelastic deformation most materials exhibit a resistance to slip called strain harden- 
ing. Hardening occurs due to interaction of dislocations with precipitates, grain boundaries, 
or other dislocations and is often a manifestation of dislocation pileups. Dislocation inter- 
action may result in internal stress, causing a resistance to further deformation. However, 
upon reversal of loading, the dislocations will propagate more easily and may cause yielding 
to occur at much lower applied stress level, resulting in the Bauschinger effect (Khan and 
Huang [1995]). 

Another important factor that affects hardening behavior is the grain size. For a coarse 
grain material, slip is restricted to the grain boundaries. However, for a fine grain material 
like Inconel 718, slip systems are also activated in the interior of the grains causing higher 
dislocation interaction in the material. This results in higher strain hardening. 
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2.2.1 Strengthening mechanisms in Inconel 718 


Generally, strengthening in metals is achieved by slowing the movement of dislocations. In 
a precipitate hardened material, the precipitates are dispersed throughout the grains. Dis- 
locations are pinned at these precipitates causing dislocation pileups. In order to overcome 
the precipitates the dislocations either climb over or shear through them, both of which 
require additional shear stress. 

Strengthening Phases 

Paulonis et al [1969] showed that the major strengthening phase in aged Inconel 718 is the 
disc-shaped 7 // (D022) precipitates rather than the spheroidal 7 ' (LI 2 ) phase found in most 
nickel superalloys. The unit cells showing ordering for the body-centered tetragonal (bet) 
7 W an d the face-centered 7 ' are shown in Figure 2.3. They observed that fully heat treated 
and deformed samples of Inconel 718, both at room and elevated temperatures, exhibited 
planar deformation banding, a slip mode characteristic of nickel-based superalloys. They 
pointed out that the alloy exhibits good tensile strength between cryogenic temperature to 
about 705° C, even with as low as 20% by weight of the strengthening 7 " phase. Excellent 
weldability was attributed to sluggish precipitation kinetics of the 7 " phase. They found 
that the rapid deterioration of the tensile and creep strengths from 650-760° C was primarily 
due to rapid coarsening of the 7 " phase and partial solution of the 7 ' phases with concurrent 
formation of stable orthorhombic N^Nb. Nevertheless, 7 " forms coherently and tends to 
remain so throughout the coarsening period. Coherency of the precipitate is a significant 
factor that contributes to the strength of the alloy. 

Deformation Mechanisms 

Analysis of deformed samples is difficult because the combination of high particle density, 
coherency strain contrast, and dislocation strain contrast results in extremely complicated 
electron micrographs (Paulonis et al [1969]). Compression at room temperature causes {1 
11} planar slip banding. Tensile loading at elevated temperature results in heterogeneous 
slip which is confined to planar bands. Precipitate shearing within the slip bands is also 
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Figure 2.3: Unit cells showing ordering: (a) bet ( 7 ") structure (b) fee ( 7 ') structure 


observed during tensile deformation. 

Oblak et al [1974] observed that since 7 " particles produce a tetragonal distortion of the 
matrix, the specific variants of 7 " can be controlled by the application of stress during 
aging. To restore order in the DO 22 phase, atoms have to displace differently in different 
directions. Also, dislocation movement to restore order in different variants of 7 " is different. 
Hence order strengthening is a complex interplay between the stress axis, Burgers vector of 
dislocations and the variant present in the microstructure. Aging the material under tensile 
and compressive loads suppresses or enhances the formation of a particular variant of D 022 - 
This results in a stronger material when aged under tensile load and a weaker material when 
aged under compressive load. This suggests that different variants of ')" interact differently 
with the applied stress field and this could be a possible cause of the SD. 

While no evidence was found for a stacking fault mode of shear, motion of dislocation pairs 
seemed to cause shear (Oblak et al [1974]). Smaller particle sizes result in dislocations being 
arranged in quadruplets rather than in pairs. While the lowest energy mode of deforma- 
tion is the motion of dislocation quadruplets, the comparison with experimental data on 
increase in critical resolved shear stress (CRSS) points towards deformation by motion of 
dislocation pairs. For these reasons, it is felt that while order ensures the pairwise motion of 
dislocations, the principal obstacles to these pairs are not the ordered particles per se, but 
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rather their coherency strain fields. Hence, a coherency strengthening seems more probable 
than order strengthening. Also, the temperature dependence of CRSS is better predicted 
by coherency strengthening model. A possible role of the 7' phase may be to inhibit de- 
formation by stacking fault mode of shear. Hence, extremely stable microstructures can be 
obtained by formation of compound 7 r/ -7 ; precipitate particles. 

Small and Large Deformation Effects 

Amongst other mechanisms in Inconel 718 , the possibility of shearing ordered precipitates 
by deformation twinning is discussed by Sundararaman et al [ 1988 ] for room temperature 
deformation. Circular dislocation loops were found near 7' and a large number of dislocation 
pairs were found near 7" phase in samples deformed to 2% strain. Some dislocations were 
observed to split up into partial dislocations on entering 7" precipitates. Offset produced 
on 7" particles due to the passage of the deformation bands indicates that these precipitate 
regions contain multi-layered stacking faults. Also observed were two sets of deformation 
bands, lying on two different {111} plane variants. 

For large strains ( 19 %), Sundararaman et al [ 1988 ] observed similar deformation bands in 
the material. In addition, the presence of deformation twins was seen within these bands. 
The flow stress required for deformation twinning is less than that required for precipitate 
bypassing and hence twinning is favored. Deformation of samples aged for longer periods 
resulting in coarser precipitates showed that the 7" particles had apparently fragmented 
into slices along { 111 } planes. Deformation twinning of 7" precipitates conserves order 
within the twin bands, hence the matrix-precipitate lattice correspondence for the twinned 
regions is the same as that for the twin-free regions. Some matrix twin segments appeared 
to propagate through a few 7" precipitates. Offsets caused by the deformation twinning 
leads to rotation of the macroscopic 7" habit plane. Each precipitate variant is associated 
with a favored twin-plane (of the type {1 1 1}), the selection of which is presumably dictated 
by the relative orientation of the stress axis. Again, deformation twinning is not seen in 
samples with very fine sized 7" precipitates. 
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Work Hardening 


The work hardening rate in precipitate hardened alloys generally increases abruptly as the 
precipitate size grows beyond a certain critical value (Sundararaman et al [1988]). This 
effect is attributed to the transition from precipitate shearing to the precipitate bypassing 
mechanism regime where the deformation results in the accumulation of dislocations in the 
vicinity of the non-deformable precipitates. In Inconel 718, however, a drop in the work 
hardening rate is observed when the precipitate radius is greater than 10 nm. This strongly 
suggests that precipitate bypassing mechanisms, like Orowan, are not active. The passage 
of shear through the precipitate particles by deformation twinning does not increase the 
dislocation density to any significant extent and consequently this deformation mechanism 
is not expected to increase the work hardening rate. The stress required to nucleate defor- 
mation twins in the precipitate particles is larger than that required for the growth of such 
twins. Thus, nucleation is immediately followed by easy propagation of deformation twins 
across the precipitates. This brings down the work hardening rate when precipitate size is 
greater than 10 nm. 

Effect of Strain Rate and Temperature 

The effect of strain rate and temperature on the thermo-mechanical behavior of Inconel 
718 was investigated by Zhou et al [1993] to understand the relationships between the dy- 
namic restoration process and the flow stress-strain curves. Also, changes of microstructure 
associated with different degrees of reduction during compression processing were investi- 
gated. Study of stress-strain response shows higher peak stresses at lower temperatures. 
The flow stress is mildly sensitive to strain rates between 0.1-0.05 sec -1 , but this sensitivity 
increases with temperature. Dynamic recovery and dynamic recrystallization are observed 
in the microstructure throughout the range of the experiments. Nucleation of dynamically 
recrystallized grains occurs at strains corresponding to peak stress (referred to as peak 
strain) and the dislocation density in the deformed region is considerable at this strain. 
The fraction of the recrystallized grains increases with temperature. Preferred nucleation 
sites for the recrystallized grains axe the strained grain boundaries. Dynamic recovery is 
a slow process and hence the dislocation density can increase considerably before dynamic 
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recrystallization initiates at peak strain. 


A critical dislocation density, which depends on the condition of deformation, must be 
exceeded before recrystallization can proceed. An increase in temperature leads to an 
increase in grain boundary mobility and a decrease in the critical dislocation density for 
dynamic recrystallization, which in turn results in a decrease in peak strain. The strain 
rate has the opposite effect. The fraction of dynamically recrystallized grains increases with 
temperature, and the fully recrystallized microstructure is obtained at temperatures greater 
than 1050°C and a strain of around 0.7. Beyond this peak strain, more dislocations are 
generated which is a driving force for further nucleation and grain growth. 

2,2,2 Mechanisms Causing Strength Differential Effect 

Several hypotheses were proposed to explain the strength differential (SD) effect. These 
include microcracking, residual stresses, internal Bauschinger effect, particle-dislocation in- 
teractions, and volume expansion during inelastic deformation. In this section a brief history 
of these investigations is presented. Comparing and correlating these with the deformation 
mechanisms in Inconel 718, an attempt can be made to qualitatively arrive at the most 
probable cause of the SD. 

Volume Expansion Hypothesis 

Leslie and Sober [1967] were among the first to observe SD in martensitic steels. They found 
that untempered carbon martensite is significantly stronger in compression than in tension. 
The SD effect appeared to increase with increasing carbon content. Radcliffe and Leslie 
[1969] were the first to attribute SD effect to volume expansion during inelastic deformation, 
which would theoretically lead to larger stress in compression than in tension. The volume 
expansion hypothesis is also supported by Drucker [1973] for plastics and granular media. 

In some cases the effect of hydrostatic stress is to cause an increase in dislocation density, 
which can result in permanent volume expansion. However, some hydrostatic pressure de- 
pendent materials do not exhibit a permanent volume change during deformation. Spitzig 
et al [1975] have reported that the yield strength of quenched AISI 4310 and 4330 steels 
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increases evenly in both tension and compression with increasing hydrostatic pressure. Fur- 
thermore, they detected a permanent volume change proportional to inelastic strain. They 
suggested that the change in volume is due to a large increase in the dislocation density. 
Consequently, SD of 6% is observed for both steels, which remained constant with increase 
in inelastic strain and hydrostatic pressure. They proposed a modified yield function to 
account for this SD (see equation 3.43). 

Precipitate Dislocation Interaction 

Rauch [1975] speculated that a significant generation of dislocations would lead to a volume 
change, especially in aged materials where pre-existing dislocations are immobilized by 
precipitates. However, their experimental results for tempered AISI 4310, 4320 and 4330 
steels disagree. No evidence of permanent volume expansion was observed. SD appeared 
to increase with decreasing temperature. 

Rauch et al [1975] point out that though the initial density decrease (volume increase) 
in as-quenched AISI 4310, 4320, and 4330 steels can be attributed to transformation of 
retained austenite, the rate of density change continues to be high even at large strains. 
Thus SD at large strains can be better explained by increase in dislocation density during 
flow. According to them, the volume expansion hypothesis is a potential explanation for SD 
effect. Pinning of transformation dislocations by carbon and rapid mobilization of newly 
created dislocations result in volume increase during deformation. This is consistent with 
no SD effect in unaged maraging steels, which have no precipitates. All this points to 
hydrostatic stress dependence of yield and flow in martensitic steels. 

Nonlinear Elastic Interaction Hypothesis 

Another theory for the SD effect was proposed by Hirth and Cohen [1970]. They ruled 
out: (1) microcracks, because ultra fine grained martensitic steels have enhanced fracture 
toughness, are resistant to microcracking, and yet exhibit SD; (2) residual stresses, be- 
cause at sufficiently high strains, residual stresses are overcome and yet the SD persists; 
(3) retained austenite; and (4) internal Bauschinger effect, because of randomness in the 
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orientation of the martensite plates. They attribute the SD effect to the solute-dislocation 
interaction. The barrier for the dislocation segments near an interstitial carbon atom is 
influenced by the normal stress across the glide plane that is increased by compression and 
decreased by tension. Severe distortion of the iron lattice around the interstitial carbon 
atom leads to local elastic strains that become non-linear. This inevitably alters the atomic 
force-displacement relationships such that they are different in tension and compression. 
This hypothesis suggests that a significant portion of the SD phenomenon may be due to 
nonlinear elastic interactions. They were able to account for a SD of up to 3-6% using this 
hypothesis in martensitic steels. 

Kalish and Cohen [1969] suggest that displacement of iron atoms are so large that associated 
elastic strains are nonlinear and the interstitial-dislocation binding energy becomes greater 
in compression than in tension. Coherency strains around the precipitated particles result 
in greater dislocation interaction under compression than under tension. This suggestion 
is supported by Chait [1973], who tested three titanium alloys in the aged condition and 
found that coherent precipitation contributes to strength differential effect. Two of these 
alloys that exhibited large SD, possessed coherent a — /3 phase, while one alloy showed very 
small SD due to a phase in it being not coherent. Chait attributed the SD effect in Zircalloy 
to formation of deformation twins, which is more difficult in compression than in tension. 
Similar results in titanium alloys axe reported by Winstone et al [1973]. The SD effect is 
much smaller in alloys in which the precipitates are not coherent with the matrix. 

Pampillo et al [1972] add to the theory of Hirth and Cohen [1970] by suggesting that the 
nonlinear elastic behavior leads to a change in elastic moduli at sufficiently large elastic 
strains. The change occurs such that the elastic modulus decreases in tension and increases 
in compression. As a result, the internal stress, which depends on the interaction between 
the long range stress fields of the dislocations and that of the precipitates, decreases in 
tension and increases in compression, and this causes a SD. This theory explains the effect 
of carbon concentration and dislocation density for martensitic steels. Also, it explains the 
increase in SD at low temperatures. It is also consistent with the observation of Rauch 
and Leslie [1972], who reported that the elastic modulus of martensitic AISI 4320 steel was 
consistently 1-3.5% less in tension than in compression. 
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Gil et al [1999a] observed an increase in instantaneous stiffness of aged Inconel 718 in com- 
pression. They speculated that this stiffening could be associated with nonlinear interaction 
between precipitates and dislocations. However, stiffening was not observed in tension. 

Effect of Interface Decohesion 

Olsen and Ansell [1969] proposed a different theory for the SD effect in two phase alloys. TD- 
Nickel with 2% by volume dispersion of spherical TI1O2 particles tested at room temperature 
showed a 0.2% offset yield strength that was 30% higher in compression than in tension. 
The dislocations bypassed the Th02 particles both in tension and compression and hence 
is not the cause of SD. They suggested that voids form at the particle matrix interface as a 
result of interface decohesion during tensile loading prior to yielding. Interface decohesion 
prevents build-up of residual dislocation loops and high shear stress on particles, permitting 
a sufficient amount of dislocations to bow between particles to give macro yielding at lower 
stress level. This is not seen in compression. Thus, the mechanism in tension is bowing 
of dislocations and that in compression is shearing of precipitate particles, for which the 
required stress is greater and hence the SD. 

This theory applies to two-phase alloys with a weak bond between the particles and the 
matrix. Olsen and Ansell [1969] also studied AI-AI2O3 which is known to exhibit very good 
bonding characteristics. Here, no decohesion was seen and this results in no SD. Contrary 
to this, Mannan and Rodriguez [1973] observed that the SD in zirconium alloy increased 
with increasing interstitial content, even though no decohesion was observed. Here, the SD 
is due to the differences in the long range internal stress. This supports our contention that 
the SD phenomenon is due to different mechanisms in different materials. 

Stress Interaction with Precipitates 

Plietsch and Ehrlich [1997] studied the strain-controlled tension/compression behavior of 
pseudoelastic shape memory alloys of Ni-Ti type. They proposed an explanation of SD 
based on the generation of stress induced martensite (SIM) during loading. For a given 
stress axis, the largest possible transformation strains in tension loading can be more than 
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twice as high as the respective compressive strains. They studied austenitic, pseudoelastic, 
and purely martensitic samples (pseudoelastic samples were obtained by special thermo- 
mechanical conditioning and are capable of completely recovering deformation strains up 
to 8% by reversible transformation). The austenitic samples showed no SD. In martensitic 
samples, martensite is present as 12 variants that are internally twinned. During tension 
or compression loading, these variants are altered by successive transformations to a single 
variant that is capable of compensating for the applied stress most effectively. The strain 
induced by martensitic transformation is much larger in tension than in compression for 
the same stress level. This results in a moderate SD. Pseudoelastic samples on the other 
hand, exhibited maximum SD. This is because the variant group formed in compression can 
deliver only a small transformation strain. This results in a much steeper rise in stress and 
results in a large SD. 


2.3 Effect on Macroscopic Response 


At this point, the strengthening mechanisms in aged Inconel 718 and other mechanisms 
that result in a SD effect in metals, are reviewed. It is required to relate these mechanisms 
to the macroscopic response (experimentally observable). 

Drucker [1973] ruled out the possibility of microcracks or residual stresses being the cause 
for SD effect. According to him, it is the inhibiting influence of compression normal to 
the slip plane or twin plane which causes SD. Thus, it is most critical to verify whether or 
not the application of hydrostatic pressure (which adds equal compression normal to each 
plane) raises the magnitude of yield or flow strength in shear, tension or compression. For 
a material exhibiting a SD, a permanent volume change is a necessary accompaniment to 
inelastic deformation if the normality flow rule is applied. Generation of dislocations as well 
as vacancies in metals is associated with a volume increase and so might lead to similar 
macroscopic behavior. Stability at a material point implies convexity of the yield surface 
and for convexity to hold, the effect of added hydrostatic pressure on solid materials is 
limited by the theoretical strength. Drucker concludes that the most critical experiment 
to understand the SD is a tension test or a compression test in the presence of hydrostatic 
pressure. 
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The pressure dependency of flow results from the fact that the basic flow events like dislo- 
cation glide are pressure dependent (Spitzig [1979], Spitzig and Richmond [1984]). Spitzig 
[1979] showed that the normal stress on the slip planes affects the dislocation interactions 
in iron single crystals. The results suggest that hydrostatic pressure retards the generation 
of mobile dislocations. Jung [1981] presented a thermodynamics based model which shows 
that pressure should increase the interaction between dislocations. He used this model 
to explain dislocation annihilation and show that the pressure induced friction impedes 
dislocation motion and completely compensates for the increased interaction between dis- 
locations. Spitzig and Richmond [1984] were able to directly relate equation 1.7 to Jung’s 
model and conclude that the pressure dependence of the flow stress does not come from the 
plastic dilatancy as required by the associated flow law. Rather, it comes from the effect of 
pressure on dislocation motion. 

Jesser and Kuhlmann-Wilsdorf [1972] found that the dislocation structure of commercially 
pure Nickel is not altered by hydrostatic pressure. They conclude that pressure causes a 
strong rise in frictional stress, which is associated with the Peierls-Nabarro stress being 
pressure dependent. They relate the pressure dependence to the extra volume that the 
atoms in the dislocation core occupy. 


2.4 Possible mechanisms causing SD in Inconel 718 

As mentioned earlier, microstructural analysis of deformed samples is difficult due to high 
particle density, coherency, and dislocation strain contrast. Thus, it is worthwhile to identify 
the possible mechanisms that cause SD in Inconel 718. 

Prom the above discussion of SD, the experimental results on different materials suggest that 
a microstructure containing precipitate particles is a prerequisite for SD effect. The exact 
mechanism responsible for a SD differs depending on the microstructure of the material. 
Theories such as the internal Bauschinger effect and microcracking hypothesis, suggest that 
the SD will increase with increasing inelastic strain, which was not observed in experiments. 

Interface decohesion between the 7 " precipitates and the matrix can be ruled out because 
the y" platelets remain coherent up to large strains and also during coarsening. Thus, this 
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is not the cause for a SD in aged Inconel 718. There is no significant difference between 
the elastic moduli in tension and compression. SD persists at strains beyond 2% where 
the inelastic strains are dominant and nonlinear elastic interactions are negligible. Thus, 
nonlinear elastic interaction is not the sole cause of SD. Gil et al [1999a], however, identified 
this mechanism as a cause for stiffening in compression, which occurs at the onset of yielding 
where the inelastic strains are low in comparison with the nonlinear elastic strains. 

Differential formation of deformation twins in tension and compression does not seem to 
be a primary cause for a SD. The effect of loading on different variants of y" precipitate, 
which have different favored twin planes, can cause a difference in response under tension 
and compression. Though no definite evidence has been seen for or against this possibility, 
deformation twinning seems to be a secondary effect and is not likely to be a dominant 
cause for SD in Inconel 718. 

Two theories that seem more convincing are: (1) volume expansion hypothesis; and (2) the 
particle-dislocation interaction hypothesis. In most cases, both of these underestimate the 
observed SD effect, suggesting that both of these mechanisms contribute to this effect. It 
remains to be seen whether hydrostatic pressure has an influence on the SD in Inconel 718. 
As suggested by Drucker, tensile and compressive experiments under varying hydrostatic 
pressure are planned on aged Inconel 718. In general, hydrostatic pressure impedes the 
motion of dislocations causing an increase in dislocation density which in turn causes a 
volume expansion. 

In the event of hydrostatic pressure-independence, the most probable cause for SD in aged 
Inconel 718 appears to the interaction between coherency strain fields around the 7 " precip- 
itates and the dislocations. It appears that the precipitate-dislocation interaction causes an 
increase in the dislocation density, the coherency strains seem to cause greater dislocation 
interaction under compression than under tension. If the material is hydrostatic pressure- 
dependent, the SD can be attributed in part to the particle-dislocation interaction and in 
part to the volume expansion. 
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Chapter 3 


Continuum Plasticity 


Starting with the definitions of stress invariants, this chapter presents the basic elements of 
rate-independent (classical) plasticity theory, viz. yield criterion, flow law, hardening law, 
and loading criteria. Then, the internal variable approach to rate-dependent plasticity is 
presented which leads into the internal variable theory of viscoplasticity. Finally, a review of 
general yield functions and flow laws that go beyond classical plasticity theory, is provided. 

Any material deforms when subjected to external forces. The deformation is elastic if it 
is reversible and time-independent, viscoelastic if it is reversible and time-dependent, and 
plastic if it is irreversible or permanent. Brittle materials like glass and concrete exhibit 
very little inelastic deformation before fracturing. On the other hand, metals can undergo 
significant inelastic deformation before failure and therefore are ductile materials. 

The theory of plasticity deals with the stress-strain and load-deformation relationships 
for a ductile material beyond the elastic limit. The establishment of these relationships 
generally follows two steps: (1) the experimental observation and (2) the mathematical 
representation. The stress states that are normally achieved in an experiment are uniform, 
but the ultimate goal of any plasticity theory is a general mathematical formulation that 
can predict the inelastic deformation of materials under complex loads. 

The theories of plasticity can be established in two categories: (1) mathematical theories and 
(2) physical theories. Mathematical theories are formulated to represent the experimental 
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observations as general mathematical formulations. These theories do not require deep 
knowledge of physics of plastic deformation and are based on hypotheses and assumptions 
from experimental results. Therefore, the mathematical theories are referred to as phe- 
nomenological theories. The physical theories on the other hand, attempt to quantify the 
plastic deformation at the microscopic level and explain why and how the plastic deforma- 
tion occurs. The movements of atoms and the deformation of single crystals and grains are 
important considerations. Here, metals are viewed as aggregates of single crystals or poly- 
crystals and their response to applied loads are derived from their building blocks, namely 
single crystals and grains (Khan and Huang [1995]). 

Most applications, such as metal forming or structural design, are based on macroscopic 
quantities. Any plasticity theory must therefore contain mathematical variables that can 
be measured during experiments at the macroscopic level. To make reasonable hypotheses 
and assumptions on the basis of experiments and to fully understand the meaning and 
limitations of a proposed theory, knowledge of material structure and physics of inelastic 
deformation is very helpful. 


3.1 Invariants of the stress tensor 


The stress state at any material point may be represented by the stress tensor Oij. Its 
components in the cartesian coordinate system can be represented by a matrix of second 
order 


G xx 

G X y 

Gxz 

Gyx 

Gyy 

°yz 

G zx 

°zy 

Gzz 


(3.1) 


Since the stress tensor is symmetric, only six stress components G X x^G X y^a xz ^ayy,a yz ,a zz 
axe independent. Thus, six independent stress components determine a stress state uniquely 
and visa versa. Using the above stress tensor, the three principal stresses can be determined 
using the characteristic equation 


Wij - crSij\ = 0 


(3.2) 
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which can be expanded to 


<7 3 -/i<7 2 + J 2 <7 -/ 3 = 0 (3.3) 

where $ij is the Kronecker delta, Ji, / 2 , and / 3 are the invariants of the stress tensor. In 
terms of the principal stresses, <7i, <72 and o%, the invariants are 

1 1 = 0[ H - <72 (J3 

h = <7l <72 + <7 2 <73 + a3<7i (3.4) 

/3 = aiders* 

In plasticity theory, it is customary to decompose the stress tensor as 

Oij — O m Sjj + Sij (3*5) 

where a m is the hydrostatic stress given by 

Ojii — “ 7^{Oxx 4" Oyy + O zz ) = ~(<7i + <7 2 + < 73 ) (3*6) 

and Sij is the deviatoric stress tensor. Using the deviatoric stress tensor, the deviatoric 
stress invariants can be defined as 

Ji = Su — 0 

h = = 1[((7! - <7 2 ) 2 + (<x 2 -<t 3 ) 2 + (<7 3 -<Tl) 2 ] (3.7) 

h = \s i:i S jk S ki = 5!5 2 5 3 . 

The three invariants, /j, J 2 and J 3 were discussed in Chapter L Invariants of the effective 
stress which are used in the proposed model are defined in a similar manner (Chapter 4). 


3.2 Rate Independent Plasticity 

In this section, the fundamental aspects of classical or rate- independent plasticity are ex- 
plained. The three basic elements of classical plasticity are yield criterion, flow rule and 
hardening rule. 
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3.2.1 Yield Criteria 


In a uniaxial stress state, the elastic limit is defined by the yield stress. Under combined 
stresses, the elastic limit is defined by a surface in stress space. Mathematically, the elastic 
limit, for a general anisotropic material is expressed as 

= 0. (3.8) 


Invariant representation of yield 

For an isotropic material, the orientation of principal axes is immaterial and hence the 
principal stresses, <7i, <72 and <73 are sufficient to uniquely describe the state of stress. The 
principal stresses form the integrity basis (Spencer and Rivlin [1962]). It is customary to 
use Ji, J2, and J3 as an integrity basis. Hence, the yield function becomes 

/(/i,J 2 ,J 3 )=0. (3.9) 

As a further refinement to this criterion, evolution of material state can be incorporated in 
the integrity basis to represent the flow behavior. 

Pressure independent yield criteria 

Since yielding of most metals is found to be insensitive to hydrostatic pressure, it follows 
that shear stresses must control the yielding of such materials. There are several shear 
stress based yield criteria. Traditionally, Tresca and von Mises criteria are widely used for 
defining yield in metals. 

1. Tresca Yield Criterion states that yielding occurs at a material point, when the 
maximum shearing stress at that point reaches a critical value k. In three dimensional 
principal stress space, the yield surface is a regular hexagonal cylinder with its axis 
along the hydrostatic stress axis. 

2. von Mises Yield Criterion states that yielding occurs when the strain energy of 
distortion in a material reaches a critical value. In three dimensional principal stress 
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space, the yield surface is an infinite circular cylinder with its axis aligned with the 
hydrostatic stress axis. 

The Tresca and Mises yield functions are independent of: Ji, implying hydrostatic pressure- 
independence; and J 3 , meaning no asymmetry between tensile and compressive stresses (no 
SD). 


Pressure-dependent yield criteria 


Geological materials such as soils, rocks and concrete are hydrostatic pressure sensitive. 
Hence, their yield representation must involve I\. Some metallic materials are also pressure 
sensitive and require a different yield representation compared to the conventional Tresca 
or Mises criterion (Drucker [1973], Spitzig [1975]). The simplest and most popular two- 
parameter models are those of Mohr-Coulomb and Drucker-Prager. 

Mohr-Coulomb Criterion 

Mohr’s criterion is based on the assumption that the maximum shear stress is the decisive 
measure of yielding. It states that yielding occurs when the radius of the largest principal 
circle touches the yield envelope, which are straight lines in the (< 7 , r) space as shown in 
Figure 3.1. Mathematically, this can be expressed as 


r\ — c — a tancj) 


(3.10) 


in which c is the cohesion and (j) is the angle of internal friction ; both are material constants 
determined by experiments. In special cases of frictionless materials, for which 0 = 0, 
equation 3.10 reduces to the Tresca criterion. In terms of the stress invariants, the criterion 
can be written as 


/(Ji, e/ 2 , 0) = \l\sin<l>+ \/j 2 [sin (0 + 4 - ~-^cos(9 + ^)sin(p] - c coscj) 

3 3 y 3 3 


— 0 (3.11) 


where 6 depends on J 2 and J 3 as 


cos (30) = 


3v/3 J z 

9 3 

Ji 


Graphically, this criterion is shown in different stress planes in Figure 3.1. 
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Figure 3.1: The Mohr-Coulomb Yield Criterion 


Drucker-Prager Criterion 

The Drucker-Prager criterion is a simple modification of the von Mises criterion to include 
the effect of hydrostatic pressure on yielding. Mathematically, 

f{Iu J 2 ) — ctli + \[Jz — fc = 0 (3.12) 

where a and k are material constants. When a is zero, it reduces to the von Mises criterion. 
The yield surface here, is a right circular cone in the principal stress space, as shown in 
Figure 3.2. 

Unlike the Mohr-Coulomb hexagonal yield surface, the Drucker-Prager surface is smooth. 
This facilitates its application in plasticity theory because it becomes easy to calculate the 
normal vector to the yield surface numerically. 
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Figure 3.2: The Drucker-Prager Yield Criterion 


3.2.2 Flow Law 


Plastic Potential Theory 


The general mathematical treatment of the constitutive equation for plastic deformation 
was proposed by von Mises in 1928. He noticed that in elasticity theory, the strain tensor 
was related to the stress through an elastic potential function, the complementary strain 
energy U c such that 


£ ij 


dU c 

dc?ij 


(3.13) 


By generalizing and applying this idea to plasticity theory, Mises proposed that there exists 
a plastic potential g{a t j ), such that the inelastic strain rate ejj could be derived from the 
following flow rule: 


e ij = A 


dg 

da. 




(3-14) 


where A is a proportional positive scalar factor which can be determined from the yield 
criterion. Plasticity theory based on the above flow rule is called plastic potential theory. 
The following remarks should be noted about the above flow rule: 


1. Geometrically, the plastic potential g{o%j) = 0, represents a surface in the stress space 
and efj can be represented by a vector in this space. The inelastic strain rate vector 
is normal to g(aij) = 0. Therefore, equation 3.14 is also referred to as the normality 
flow rule in plasticity theory. 
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2. If inelastic flow is hydrostatic pressure independent, and thus incompressible, the 
corresponding plastic potential surface in the principal stress space must be cylindrical 
(not necessarily with a circular cross section) with a\ — 02 — 03 as its axis. ejj is 
perpendicular to this axis representing only shape changes of the surface. This is not 
the case for hydrostatic pressure-dependent materials. 

3. A common approach in plasticity theory is to assume that the plastic potential func- 
tion g{oij) is the same as the yield function so that the flow law can be written 

as 

4 “ ^ (315) 

and the inelastic strain rate vector is normal to the yield surface. This is called the 
associated flow rule. On the other hand, if g / /, the flow rule is called nonassociated. 

4. In general, experimental observations show that inelastic deformation of metals can 
be characterized quite well by an associated flow rule, but for some porous materials 
a nonassociated flow rule provides a better representation of inelastic deformation. 

Prandtl-Reuss Incremental Stress-Strain Relations 

The Prandtl-Reuss flow theory is based on the following three assumptions: 

1. The principal axes of the inelastic strain increment are coincident with those of the 
current stresses. 

2. The inelastic strain increments are proportional to the deviatoric stress tensor. 

3. No volume change occurs during inelastic deformation. 

4. Initial response is isotropic. 

In contrast to the deformation theory, the stress-strain relation in the flow theory is given 

in an incremental form. Here, the inelastic (deviatoric) strain increments are given by 

de\j = dX Sij (3.16) 
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Figure 3.3: Graphical Representation of Flow Rule 

where dX is a positive scalar factor. This relation is called the Prandtl-Reuss relation. An 
interesting observation is that the associated flow law for the Mises yield function is the 
Prandtl-Reuss relation. This means that the Prandtl-Reuss relations are associated with 
the Mises yield criterion. The flow rule is graphically shown in Figure 3.3. 

Subsequent Yield Surface 

For an elastic-perfectly plastic material, the subsequent yield surface is the same as the 
initial yield surface since there is no increase in stress carrying capacity of the material 
beyond the yield stress. But, for hardening materials, the two surfaces are different and 
there are various models to describe how the initial surface evolves into the subsequent yield 
surface. Before discussing these models, the criteria for loading and unloading should be 
defined. 

Loading criteria 

For a hardening material, if the stress state tends to move out of the yield surface, we have 
a loading process and elastic-plastic deformation is observed. Additional plastic strains are 
produced and the configuration of the current yield or loading surface changes, so that the 
stress state always stays on the subsequent loading surface. If the stress state tends to move 
into the yield surface, we have an unloading process. Only elastic deformation occurs and 
the loading surface remains unchanged. The other possibility of stress change from a plastic 
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e 

Figure 3.4: Loading Criterion for a Work-Hardening Material 


state is that the stress point moves along the current yield surface. This process is called 
neutral loading. Mathematical expressions for classifying these criteria are called loading 
criteria , and may be expressed by: 

/ = 0 , don > 0 : Loading 

d(Tij 

f = 0 , -^—doij = 0 : Neutral loading (3.17) 

doij 

r\ ft 

f = 0 , — — daij < 0 : Unloading. 

OOij 

The gradient vector, nij{=d f / doij) is the outward normal to the yield surface / = 0. These 
criteria are illustrated schematically in Figure 3.4. 


Consistency Condition 

During inelastic deformation, the elastic region has to change so that the stress state remains 
on the yield surface. Mathematically, this condition is expressed as 

f(o + da, k + dn) = 0 (3.18) 

where n is the hardening parameter. In incremental form, it may be rewritten as 

% da + - 0. (3.19) 

do OK 

This equation is known as the consistency condition , and lends itself to the determination 
of the magnitude of plastic strain increment. The consistency condition ensures that the 
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plastic loading path begins at the current yield surface and ends at the subsequent surface, 
and that changes in the size and location of the yield surface are consistent with the adopted 
hardening rule. 

Hardening Laws in Plasticity 

For work hardening materials, stress states can exist beyond the initial yield surface and the 
yield surface evolves such that the current stress state always lies on it. There are various 
models describing this evolution of the yield surface. These are called hardening models. 

Isotropic Hardening 

This is the simplest hardening model and is based on the assumption that the yield surface 
expands uniformly without distortion, as the inelastic flow occurs, as shown in Figure 3.5(a). 
If the yield surface has the form /(o^j) = /c 2 (e p ), the size of the yield surface is governed 
by the value of A: 2 , which depends upon either effective inelastic strain or plastic work. 
Isotropic hardening can be observed in materials with forest dislocations. 

Since the loading surface expands uniformly (or isotropically), it can not account for the 
Bauschinger effect exhibited by most structural metals. Hence, it will not lead to realistic 
results when complex loading paths involving significant changes in the direction of stress 
vector are considered. 

Kinematic Hardening 

The kinematic hardening model assumes that during inelastic deformation the yield surface 
translates as a rigid body in the stress space, keeping its shape, size, and orientation the 
same as that of the initial yield surface. Kinematic hardening is often associated with 
dislocation pileups in a material. This model has the form 

/(oy-oy)-* 2 = 0 (3.20) 

where A; is a constant and otij axe the coordinates of the center of the yield surface, which 
changes as inelastic deformation progresses. This is shown graphically in Figure 3.5(b). 
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Figure 3.5: Hardening Rules: (a)Isotropic; (b)Kinematic 


Mixed Hardening 

A combination of isotropic and kinematic models leads to a more realistic mixed hardening 
model whose loading function can be expressed as 

f(o ij -(x ij )-k 2 (e p )= 0. (3.21) 

In engineering application for metals, the concept of mixed hardening is attractive. The 
loading surface translates and uniformly expands in all directions; i.e. it retains its original 
shape. With mixed hardening, different levels of Bauschinger effect can be simulated. 

Drucker’s Stability Postulate 

Drucker [1959] defined a stable material as one which satisfies the following conditions: 

1 . The work done by an external agency during the application of an added set of forces 
on the change in displacement it produces, is positive. 

2. The net work performed by the external agency over the cycle of application and 
removal of the added set of forces on the changes in displacements it produces is 
nonnegative. 
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These criteria are called Drucker's Stability postulate . The first is for a work increment 
during loading and is called stability in small and the second is for a work increment after 
full load/unload cycle, and is termed stability in cycle. Mathematically, these reduce to 

dade > 0 loading ; (3.22) 

dode 1 > 0 complete cycle. (3.23) 

Stability in a complete cycle leads to the following: 

(<7ij - o*j)e!j > 0 (3.24) 

where Oij is a stress state on the yield surface and o* 3 is any possible stress state inside or 
on the yield surface. This is called maximum plastic dissipation postulate. 

Convexity and Normality 

If the plastic strain coordinates are superimposed upon the stress coordinates, as shown in 
Figure 3.6, equation 3.24 can be interpreted geometrically as the scalar product of the stress 
increment vector [o%j — o* 3 ) with the strain increment vector de 1 - . A positive scalar product 
indicates an acute angle between these two vectors. The stability postulate therefore leads 



Figure 3.6: Stability in Cycle: Stress Path ABC produced by External Agency 


to the following consequences for work hardening materials (Drucker [I960]): 

1. Convexity: The initial yield surface and all the subsequent yield surfaces must be 
convex. 
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2. Normality: The inelastic strain increment vector, dejj must be normal to the yield 
surface at a smooth point and must lie between adjacent normals at a corner. 


3.3 Rate Dependent Plasticity 

One way that loading history affects the constitutive relation is through rate sensitivity ; 
the deformation produced by a slow loading rate is different, almost invariably greater than 
that produced by rapid stressing. A particular manifestation of the time dependence is the 
fact that the deformation will in general, increase in time at a constant stress; called creep. 
Relaxation , on the other hand is the decrease in stress with time for a given strain. 

Many physical and chemical processes in metals are thermally activated and hence the 
dependence of material behavior on temperature. For most metals, creep is thermally 
activated if the temperature is higher than 30% of the melting temperature. Creep can 
cause failure, either due to excessive deformation or creep rupture. These processes are 
often governed by the Arrhenius rate equation , which has a general form 

e ss = Ae^ kT (3.25) 

where e ss is the steady state creep rate, A E is the activation energy, k is Boltzmann’s con- 
stant and T is the absolute temperature. The rate sensitivity of a work- hardening material 
itself increases with temperature. The above equation permits, in principle, the simul- 
taneous representation of the rate sensitivity and the temperature sensitivity of metallic 
materials. 

An alternative way of representing the dependence of material response on temperature and 
rate is by assuming that the strain, in addition to stress and temperature, depends on an 
array of variables, ap. These variables are called internal , or hidden variables, that usually 
take on scalar or second rank tensor values. The strain is accordingly given by 

e = e(<7 ij ,T,a /3 ). (3.26) 

Another important effect observed in metals is creep which depends on stress level and 
temperature. Rate independent plasticity does not represent these effects. This section 
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starts with general principles in mechanics that lead to the internal variable theory of 
plasticity. 

3.3.1 General Principles in Mechanics 

Before we discuss the internal variables approach, it is necessary to lay out the universal 
laws that must be obeyed during any process, regardless of the properties of the material 
which is undergoing the process. There are five such principles to be considered in plasticity 
theory (Khan and Huang [1995]). 

Conservation of Mass 

According to the law of conservation of mass , no mass can be created or destroyed in a 
given volume v of the material, which leads to 

p dv = constant (3.27) 

where p is the material density. 

Conservation of Momentum 

The rate of change of total momentum of any given set of particles equals the vectorial sum 
of all the forces acting on this set of particles. For a set of particles that currently occupies 
a spatial volume, v, with surface traction and the body forces acting on it, this law leads to 

+ pbi = pa t (3.28) 

where 6j’s are the body forces per unit volume and a t ’s are the components of the acceler- 
ation of the body. Equations 3.28 are also known as the local equations of equilibrium. 

Conservation of Angular Momentum 

This principle is also a generalization of Newton’s second law of motion for continuum 
mechanics. It states that the rate of change of moment of momentum for any given set 
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of particles equals the vectorial sum of the moments acting on it. Mathematically, the 
principle of conservation of angular momentum can be stated as 

I ip s ijk f j o,k)dv — I (ej jk r j tk)ds + I ip ^ijk ^'j bk)dv (3.29) 

J v Jr Jv 

where e^jt is the alternator, t x is the traction force vector and r t is the position vector of 
the body (volume v enclosed by surface R) under consideration. This, along with equation 
3.28, establishes the symmetry of the stress tensor, i.e. ct x] = oj x . 

Conservation of Energy: The First Law of Thermodynamics 

Extensive experimental observations indicate that energy can never be created or destroyed 
in the universe, but can only be transformed from one form to another. For a closed system 
the total rate of work done on the system by all the external agencies must equal the rate 
of increase of the total energy of the system. This principle is also called the first law of 
thermodynamics. If we consider only mechanical and thermal energy in a closed system, 
the first law can be written as 

pu = Oijiij + pr — qij (3.30) 

where u is the internal energy per unit volume, r is the internal heat source per unit 
volume, q iti is the heat flux out of the system through its boundaries, and represents 

the mechanical work done by the external forces that is not converted into kinetic energy. 

Clausius-Duhem Inequality: The second Law of Thermodynamics 

The second law of thermodynamics limits the direction of the energy transformation. When- 
ever a transformation occurs, the energies involved must obey the first law. It has been 
experimentally observed that while some energies transform from one type to another, there 
are other types of transformations that are impossible. For example, heat flow can occur 
from a warmer system to a colder system, but the reverse heat flow can never occur. The 
kinetic energy of a moving body can be converted into heat by friction, but the heat caused 
by friction can never be converted into kinetic energy. While the first law cannot describe 
these observations, it is the second law that governs this directional phenomenon observed 
in the energy transformation processes. 
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The second law postulates that there exists a state function called the entropy; the change 
in entropy is given by 

AS - ft-ft = 

AS = ft -ft > 

where dQ is the heat input during the process, T the temperature, and indices 1 and 2 
denote the starting and ending points of the process. It can be seen that the change in 
entropy of a system can never be negative; it is zero for a reversible process and positive for 
an irreversible process. If only mechanical and thermal energies are considered, the second 
law can be written as 

pi - 4 + (§). > 0 (3.31) 

where s is the specific entropy. This relation is called Clausius- Duhem inequality. Using 

equation 3.30, r can be eliminated and we get 

pTs - pit + aijtlj - —qiTi > 0. (3.32) 

By introducing a thermodynamic potential H , called the Helmholtz free energy , given by 

H = u — sT, we can rewrite equation 3.32 as 

-pH - psT + ^ qi Ti > 0. (3.33) 

In thermodynamics, the internal energy u, entropy s, heat flux q and the stress a are 
all considered state functions that can be determined by the state variables using state 
equations or constitutive equations . 

3.3.2 Internal Variables: General Theory 

In thermoelasticity, internal energy it, entropy s, and heat flux q, can be fully described by 
the current value of stress and temperature. The situation becomes more complex when 
deformation is inelastic. Here, in addition to the current values of the state variables, the 
history of the deformation is also important. This requires identifying more state variables 
to describe the deformation history. It is very hard to enumerate all of the relevant state 


for a reversible process 
for an irreversible process 
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variables from macroscopic observations alone. Microscopic deformation considerations are 
needed and some assumptions have to be made to use certain macroscopic (observable) 
variables as representatives of the microscopic phenomenon. After the state variables are 
chosen, the mathematical forms of the constitutive equations should be determined. This 
involves experimental evaluation and mathematical formalization, which at times are very 
complex (Lubliner [1990]). 

In thermomechanics, the state variables are identified by using the concept of internal 
variables. It is postulated that the current state of an inelastically deformed solid can 
be determined by the current values of stress, temperature, and a set of internal variables. 
The history of deformation is indirectly included in the evolution of these internal variables. 
Mathematically, this can be stated as 

u = u(<7ij,T, ap) 

s - s{crij,T,ap) (3.34) 

Q = q(°ij,T,ap) 

where ap, i=l,2,...,n are internal variables which can be scalars, vectors or tensors. These 
variables can be considered in two categories; (1) physical variables describing the aspects 
of local physico-chemical structure which may change spontaneously, and (2) mathemat- 
ical constructs or phenomenological variables , in which case the functional dependence of 
stress(or strain) on the internal variables, and their rate equations, is assumed a-priori. 

Thermomechanical Aspects 

An equilibrium state of a system is a state that has no tendency to change without a change 
in external controls. The local state (cr, T, a) may be called a local equilibrium state if 
the internal variables remain constant at a constant stress and temperature. This can be 
mathematically stated as 

a = gi(o,T,a) — 0, i— l,2,...n. (3.35) 

In an elastic continuum, every local state is an equilibrium state, though the continuum 
need not be globally in equilibrium. On the other hand, existence of non- equilibrium states 
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is an essential feature of rate- dependent inelastic continua. Such states evolve in time 
by means of irreversible processes, of which creep and relaxation are examples. Hence, 
the thermomechanics of inelastic continua belongs to the domain of thermodynamics of 
irreversible processes (Lubliner [1990]). 


3,3.3 Generalized Flow Potential 


A general concept of the flow potential is due to Rice [1970,1971]. First, a generalized flow 
potential is defined which is assumed to depend on the the stress, temperature and a set of 
internal variables. The flow and evolution laws can be obtained by simply differentiating 
this potential. 


For a strain formulation, generally the Helmholtz free energy, H(eij, T, e\^ A^) is used as the 
generalized potential. This potential is convex and the mechanical and thermal quantities 
can be obtained as follows: 


a ij 



dH _ dH 

df' a P~ ~ 


(3.36) 


For a stress formulation, Gibb’s complementary free energy (3>), given by 


$(aij,T,ap) = aijtij - H 


(3.37) 


is used as the generalized potential, from which the state variables can be obtained as 

d<Z> 


e ij 


do { j ’ 


d<f> , d<f> 

S = - — , A/3 = - 


dT' dap 

Further, the evolution is governed by the dissipation potential , Q(oij,T,ap) as 

A/3 = 


(3.38) 


./ on 

€ '> = 


dn 


13 


dap 


(3.39) 


3.3.4 Internal Variable Theory of Viscoplasticity 


Most researchers use the term viscoplasticity in the classical sense, that is, to denote the 
description of rate-dependent behavior with a well defined yield criterion, but this usage is 
not universal. Others, following Bodner [1968], use the term for highly non-linear behavior, 
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without a well defined yield, that is characteristic of metals, especially at higher tempera- 
tures. These are unified models. Nevertheless, both of these are subclasses of the internal 
variable models (Lubliner [1990]). 

Classical Viscoplasticity 

Suppose we have a continuous function f(a,T,a) such that there exists a region in stress 
space where /(cr, T, a) < 0 and efj = 0, then f(a, T, a) = 0 is the threshold surface in 
the stress space and the elastic region forms its interior. This definition does not require 
simultaneous vanishing of all the internal variable rates {dtp) in the elastic region. Hence, 
phenomenon like strain-aging, which require evolution of local structure while the material 
is stress-free, can be represented. However, this is of importance only for processes whose 
t im e scale is comparable to the relaxation time for strain-aging. 

The dependence of the threshold function on the internal variables ap describes the hard- 
ening properties of the material. During hardening, the threshold function decreases from 
a positive value toward zero at a constant stress and temperature, that is / < 0. Similarly, 
softening is characterized by / > 0. In the limiting case / = 0, i.e. / is independent of ap 
and the material is perfectly plastic material. 

Unified Viscoplasticity 

According to Bodner [1968], yielding is not a separate or independent criterion but is a 
consequence of a general constitutive law of the material behavior. Since the 1970’s, several 
constitutive models for the rate-dependent inelastic behavior of metals have been formulated 
without a formal hypothesis of threshold surface. These models can successfully represent 
creep especially under high temperature without a decomposition of strain into plastic and 
creep strains. They have consequently come to be known as unified viscoplasticity models, 
and are particularly useful for the description of bodies undergoing significant temperat ure 
changes. 

Perhaps the simplest unified model is due to Bodner and Partom [1972,1975] in which the 
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flow equations are given by 


ejj = 4>S tJ (3.40) 

where 

<t> = 4>(D 0 ,n,Z) (3.41) 

and D 0 , n are material parameters and Z depends on the inelastic power {W{). The rate 
equation is 

Wi = (Ti/ij = 2J 2 <t>(Wi,J 2 ). (3.42) 

Here, temperature dependence is achieved through temperature-dependent material param- 
eters. 

More recently, sophisticated unified viscoplastic models, that describe many features of the 
behavior of metals at elevated temperatures have been developed. Arnold and Saleeb [1994] 
developed a Generalized Viscoplastic model with Potential Structure (GVIPS). In GVIPS, 
the total strain is decomposed into elastic and inelastic strains. However, it is still unified in 
the sense that the inelastic strain in not further decomposed into plastic and creep strains. 
It is a fully associative, multiaxial, non-isothermal, non-linear kinematic hardening model 
that accounts for most aspects of temperature and rate-dependent inelastic deformation. 
This is the baseline model that is modified in the current work and will be presented in 
detail in Chapter 4. 


3.4 Review of Related Work 

It is clear that the SD cannot be represented by the classical Mises yield criterion with 
an associated flow rule. This prompts us to look at more general constitutive models that 
can better represent the yield and flow behavior of the materials that exhibit a SD effect. 
Reviewed here are some forms of threshold functions, definitions of flow, and flow rules that 
have been used to capture different aspects of flow behavior in materials. 
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3.4.1 General Threshold Functions 


One of the first general yield functions for metals was proposed by Drucker [1949] as de- 
scribed in equation 1.5. A function based on J 2 -J 3 was better than either Tresca or von 
Mises functions in correlating yield surface data for the aluminum alloy. Further, he proved 
that the deformation theory of plasticity falls short in representing a general state of stress 
in a material and is not compatible with the J 2 -J 3 representation of yield. 

Spitzig et al [1975] proposed a generalized yield function to account for the observed SD 
effect in tempered AISI 4310 and 4330 steels. 

/ = J 2 1/2 -a h-b J 3 1/3 - c (3.43) 

where a, 6 , and c are material constants. To be able to correctly predict the observed 
volume expansion during inelastic deformation, they suggested the use of a non-associated 
flow law. 

Lee and Ghosh [1996] addressed the problem of expressing the non-coaxiality of stress path 
in constitutive modeling. Due to shear banding, all the stress increments in a given stress 
path do not have the same principal axes. This is referred to as non-coaxiality and is 
attributed to the rigid body rotations associated with pure shear, compressibility due to 
Poisson effect and plastic dilation. To account for non-coaxiality and shear banding com- 
monly observed in deforming metals, they proposed the inclusion of J 3 in the constitutive 
model. To this end, they proposed two modifications to the Drucker-Prager criterion; one 
for materials weak in shear and another for those weak in tension. 

Lissenden et al [1999] developed a method for determining rate dependent flow surfaces for 
Inconel 718. Inadequacy of J 2 based models to represent the asymmetry in yield (between 
tension and compression) led them to use the following forms for the yield functions: 

/ = a h + b J 2 1/2 - 1 (3.44) 

/ = 6 3 J 2 /2 +cJ 3 -l 

where a, b, and c are material constants that fit to the experimental data. These functional 
for ms fit the yield surface data very well and the outward normals to the surfaces were 
found to be consistent with the experimentally determined directions of the inelastic strain 
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rate vectors. 


In order to characterize the viscoplastic behavior of geological materials, Desai and Zhang 
[1987] proposed a generalized yield function of the form 

/ = h ~ (-a I\ + 7 I\)\ ex P{P h)~ P S r ] m (3.45) 

where a, 7 , (3, n, and m are response functions and S r = /J^ 2 ■ This yield function is 

continuous in stress space. It permits hierarchical development to incorporate progressive 
complexities such as associative and non-associative responses, anisotropic hardening, strain 
softening and fluid pressure. 


Mandl and Luque [1970] analyzed a fully developed shear flow of frictional granular material. 
They justified the use of Mohr-Coulomb criterion for planar elements that have normals in 
the flow plane to describe yield in this situation. They rewrote the criterion in cartesian 
coordinates as 


+ < 7 > 


yy 


sin p + c cos p — 


— c 7 . 


yy 


\2 , _2 

) "P <Txy 


lV2 


(3.46) 


Functions based on Internal Variables: Threshold Functions 


Wegener and Schlegel [1996] studied the suitability of different yield functions for approx- 
imation of subsequent yield surfaces to capture various effects of distortional hardening. 
They chose a yield function of the form 

f = f(r,z k ,k = 1...) (3.47) 

where r is the effective shear stress, Zk are the internal variables that are tensorial in nature. 
They compared four experimental yield surfaces corresponding to different load paths from 
the work of Phillips and Tang [1972] and found that a sixth order tensor for z * gave a good 
representation of the subsequent surfaces for all the cases. 


Chaboche [1977] used inelastic strain(e 7 ) and strain like symmetric second-rank tensor A as 
internal variables. The thermodynamic force conjugates to these variables are the stress-like 
variables, R and respectively, and the assumed threshold surface is 


f(a,a,R) 



(3.48) 
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where 


J 2 — 2 ( s «j a ij ) (3.49) 

a,j being the deviator of a tJ and A : 0 is the yield strength in shear. The yield surface is 
of Mises type, capable not only of expansion (as measured by R) but also of translation 
represented by a t j which marks the center of elastic region. The invariants used in the yield 
function are effective deviatoric stress invariants. 


Robinson and Ellis [1986] assumed the following form for the dissipation potential: 


fl = k [T f F n dF + ^ f G' m dG' 

LW hJ 


(3.50) 


where /i, R , H ) n, m, k are material constants. The threshold function F depends on J 2 
and J 3 and has its origin in Drucker’s form (equation 1.5). 


F(J 2 ,J 3 )= { -^ *>* -1 


(3.51) 


and the function G ’ depends on the corresponding invariants of the internal stress and 

J3) ^ 




(jf-cjf ) 1 / 3 

K 2 


(3.52) 


The effective and internal stress invariants are defined in Chapter 4. The threshold function 
reduces to the Mises criterion for c=0, and Drucker’s model for c=-1.75. 


Hopkins [1990] developed the flow and the evolutionary laws for the above model. He 
arrived at the limits on the value of c using material stability criterion. He also determined 
functional forms for F and G’ for non- monotonic loading conditions like stress reversals. 


Janosik and Duffy [1997] pointed out that ceramic materials exhibit complex rate dependent 
thermo-mechanical behavior. To account for the SD effect and the sensitivity of these 
materials to hydrostatic pressure, they introduced 7i and J 3 in the definitions of F and G ’ 
in the Wiliam and Warnke model [1975]. These functions were in turn used in Robinson’s 
model to develop the constitutive laws for ceramics. 
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3.4.2 Flow Definitions 


At elevated temperatures, metallic alloys exhibit strong time-dependent behavior and the 
concept of yield surface, in the classical sense, generally breaks down as we move into the 
realm of viscoplasticity. However, analogous geometrically and thermodynamically based 
concepts such as surfaces of constant strain rate (SCISRs) and surfaces of constant dissipa- 
tions (SCDs) have been postulated to play the same central role in viscoplastic constitutive 
theories as yield surfaces do in classical plasticity. Physical quantities used to define the 
amount of inelastic flow in a material are defined mathematically as follows: 

1. Inelastic power is the product of the stress and the inelastic strain rate. 

2. Dissipation represents the rate of work that cannot be recovered in a process and 
mathematically is the difference between inelastic power and the product of internal 
stress and internal strain rate. 

3. Equivalent inelastic strain rate is the square root of the self product of strain rate. 

Battiste and Ball [1986] used SCISRs and Clinard and Lacombe [1988] used SCDs to describe 
the dissipation potentials at elevated temperature for monolithic materials. Lissenden et 
al [1997a] pointed out that SCISRs and SCDs need not represent the same surfaces under 
all conditions and hence needs careful consideration. A good example of this would be the 
dissipation potential in equation 3.50. For this potential, the two surfaces are different and 
become identical only if the J3 dependence is removed. Of all the surface definitions that 
represent inelastic flow, the SCDs are most meaningful from the thermodynamics viewpoint. 

An important issue is that the SCDs are not experimentally measurable. At best, ex- 
periments can measure the surfaces of constant inelastic power (SCIPs) as explained by 
Lissenden et al [2000], Hence, it is important to correlate SCIPs that can be experimentally 
determined, to SCDs that are theoretically meaningful. This subject was addressed by Iyer 
et al [2000a]. 
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3.4.3 Inelastic Behavior 


In the inelastic region, the flow law relates the incremental stress to incremental inelastic 
strain. The associated flow law, which is widely used in classical plasticity theory, was 
described earlier. Here, the focus is on the non- associated flow laws that are required for 
describing inelastic flow in some materials. 

Mroz [1963] pointed out that the vector of incremental plastic strain can rotate significantly 
for small inhomogeneity, even for a smooth yield surface, thereby leading to a non-associated 
flow rule. This inhomogeneity is caused by plastic deformation at the microscale and is a 
function of grain size, temperature, and elastic and plastic properties of the grain. Accord- 
ing to him, the separation of total strain into elastic and plastic parts is difficult due to 
inhomogeneity, which affects hardening characteristics and laws of plastic deformation. 

Spitzig et al [1975] used the yield function in equation 1.7 for describing flow in marten- 
sitic steels. An associated flow law overpredicted the volume expansion by a factor of 15. 
They proposed the use of non-associated flow law for accurate prediction of volume expan- 
sion during inelastic deformation. In Robinson’s model the flow law is associative but the 
evolution law is non-associative. 

Frantziskonis et al [1986] attributed the non-coaxiality of the plastic strain increment and 
the stress increment in granular materials to intergranular friction and changes in the phys- 
ical state of the material during deformation. They proposed a correction to the yield 
function in equation 3.45 based on the deviation of the plastic strain increment from nor- 
mality. 

F = f + h(J»£) (3.53) 

where h is the correction function that depends on the three invariants J*, i=l to 3, and 
£, which is a hardening parameter. The use of F for the yield function leads to a non- 
associated flow law, which is able to capture the effects such as volume change and stress 
path dependence. 

Mandl and Luque [1970] found that the internal kinematical constraints cause the direction 
of the principal macrostress to deviate from the direction of strain increment, though the 
granular material is assumed to be completely isotropic. This non- coaxiality is perfectly 
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compatible with material isotropy and has only minor consequences on the mathematics of 
classical plasticity. They found the angle of deviation to be the half angle of internal friction 
for a critical void ratio. Drucker’s stability postulate is not satisfied in the overall sense. 
They decompose the stress tensor into two parts; one capable of performing inelastic work 
and the other not capable of performing work. The inelastic strain increment is coaxial 
with the working part of stress, thereby indicating that a weakened form of the stability 
postulate holds. 

Another instance of application of non-associated flow laws is by Kang and Wiliam [1999] 
who used a yield function that depends on all the three stress invariants and is continuous 
in stress space except at the point of equitriaxial tension. Application to concrete resulted 
in satisfactory prediction of brittle failure mode in tension and compression failure in pure 
shear with volumetric constraint. 
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Chapter 4 


Proposed Viscoplastic Model 


In this chapter, a viscoplastic model capable of predicting a SD is proposed and subse- 
quently specialized for aged IN718 at 650°(7. This model is an extension of the Generalized 
Viscoplastic Model with Potential Structure (GVIPS), developed by Arnold and Saleeb 
[1994]. First, a background relating to the internal variables used in GVIPS is presented. 
This is followed by the theoretical framework of GVIPS and constitutive equations. The 
proposed modifications to GVIPS to account for the observed SD are presented and the 
resulting formulation is derived, giving the flow and evolutionary laws. Finally, a detailed 
procedure to estimate the material parameters in the model is proposed. 


4.1 Background 

Inelasticity exhibited by the thermomechanical response of engineering materials is related 
to irreversible thermodynamic processes . These involve energy dissipation and material 
stiffness variations due to physical changes in the microstructure. Consequently, thermo- 
dynamic arguments have often been utilized in the internal variable approach in the for- 
mulation of phenomenological constitutive laws (Coleman and Gurtin [1967]; Rice [1971]; 
Lubliner [1972,1973]; Lemaitre and Chaboche [1990]). 

Thermodynamic admissibility restrictions associated with various dissipative mechanisms 
underlying the above models reduce to the well-known Clausius-Duhem or dissipation 
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inequality. The thermodynamically based constitutive equations in these models are dis- 
cussed with respect to the above restriction. Mathematical constructs of flow or dissipation 
potential and the associated normality conditions are introduced for convenience in satis- 
fying the above constraints based on simple properties of non-negativeness and convexity 
of these functions (Ponter and Leckie [1976]; Rice [1970]; Ponter [1976], Onat and Leckie 
[1988]; Robinson and Duffy [1990]; Lemaitre and Chaboche [1990]; Freed et al [1991]). 
Such forms however, do not automatically imply the existence of total or integrated forms 
of thermodynamic potentials; for example, the Helmholtz free energy or the Gibb’s com- 
plementary free energy. If the latter is assumed a-priori, the corresponding formulation 
is termed a complete potential-based structure; on the other hand, those derived from an 
assumed dissipation potential form are referred to as incomplete potential-based models. 

The complete potential-based class of inelastic constitutive equations possesses a number 
of distinct and important attributes: 

1. They constitute the cornerstone of numerous regularity properties and bounding the- 
orems in plasticity and viscoplasticity (Ponter [1976,1979,1980]). 

2. They result in sufficiently general variational structure, whose properties can be ex- 
ploited to derive a number of useful material conservation laws (Eshelby [1951,1956]; 
Rice [1968]). 

3. The discrete form of the assumed Gibb’s potential is numerically advantageous in the 
development of efficient algorithms for finite element implementation (Saleeb et al 
[1990]; Saleeb and Wilt [1993]; Maier and Novati [1990]). 


A number of variables have been used to describe the evolution of the internal structure, 
as inelastic deformation occurs. The ones usually used in sophisticated unified theories are 
defined here. 

1. Internal Stress represents the inner stress field associated with immobile dislocations. 
It is smaller than the applied stress. 

2. Effective Stress is the difference between applied and the internal stresses. 
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3. Drag Stress is based on the dislocation motion. Increase in drag stress decreases the 
inelastic strain rate and thus increases the size of the elastic region, similar to isotropic 
hardening. 

GVIPS uses two variables; the drag stress which is non-evolving, and the internal 
stress, which is tensor ial and evolves with inelastic deformation. 


4.2 Theoretical Framework 


GVIPS applies to an initially isotropic material. It is limited to small deformations and the 
initial state of the material is assumed to be stress-free. Gibb’s thermodynamic potential 
function (<!>) is assumed to depend on stress, temperature, and an array of internal variables 
(Coleman et al [1967]; Lubliner [1972,1973])* 


In its differential form, the Gibb’s potential, $(aij,ap,T), is written as follows (Lubliner 
[1973]; Ponter [1979]): 


d$ = — eijdcrij — SdT — A pdctp 


(4.1) 


where S denotes the entropy, olq the internal state variables, Ap the thermodynamic affinities 
corresponding to a^, €ij the total strain and Oij the stress tensor. It follows from equation 
4.1 that 




* 3 


S 

A* 


daij 

cM* 

dT 

dap 


(4.2) 


are defined as the equations of state (Malvern [1969]; Lemaitre and Chaboche [1990]) and 
ap, and T are force-like thermodynamic state variables while e l} , Aj, and S are the 
corresponding conjugate displacement- like variables. The most general expression for the 
total strain rate is obtained by differentiating the first of the above equations as 


tij 


_ _ d 2 $ 

dt do ij doado r 




d 2 <f> 

doijdap 




S 2 # 

fo~dT 


T. 


(4.3) 


Two options are available for describing the flow and evolutionary equations. The first 
option assumes a fully coupled form, i.e., the inelastic strain rate is intimately linked to 
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the thermodynamic internal state, hence its functional dependence is completely defined 
once $ is assumed. This form is too restrictive, in that it requires proportionality between 
the inelastic strain and the rate of change of internal state (Arnold and Saleeb [1994]). As 
a result, the classical definition of steady state creep which requires an evolving inelastic 
strain at a constant internal state, cannot be attained using this coupled form. This leads us 
to the second option; a decoupled Gibb’s form, in which the selected internal state variables 
are grouped a priori by separation of the inelastic strain as an independent state parameter 
and suppressing all stress dependency of the remaining associated internal state groups in 
the selected $ function. This form indeed allows the classical steady state creep. 


4,2.1 Decoupled form 


In the decoupled form, the evolution of inelastic strain is independent of the internal force- 
like state variables, associated with material inelasticity. It is this separation that 
allows for inclusion of the classical notation of steady state creep. The Gibb's potential for 
an isothermal case is expressed as 

$ = E{aij) - (Tije-j + H'(ap) (4.4) 

where E is negative of the elastic strain energy, the second term is inelastic work, and H' 
is a material function. 


Differentiating equation 4.4 with respect to stress, gives the total strain rate as a sum of 
two components; elastic (or reversible) (e® ; ), and inelastic (e^). 


tij — 


d(T{j 


% 


4- e 


IJ 


(4.5) 


where 




do ij do fs _ 


and inelastic strain rate, 6^, is defined separately. 


(4.6) 


For an isothermal process, the energy balance requires that the external work rate be equal 
to the sum of the internal work rate and the dissipation, i.e. 


a ij e ij 


+ <xp A/?. 
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The second law states that the change in entropy is positive. For an isothermal process 


Q{crij,ap) - (Tijefj - agAg > 0. 
This is called the dissipation inequality . 


(4.7) 


The inelastic strain rate ejj is defined such that the dissipation inequality is satisfied. Dif- 
ferentiating the dissipation potential with respect to <t,j and ag we get the flow and 
evolution laws, 


,i dfl 
tl} = d^ 


*3 


(4.8) 


A r = 


m 

dag 


respectively. Using the Gibb’s potential, the evolution law can also be written as 


Ait = 


dt 


d<& 


da 


P\ 


Qgi dq 


where 


Qgi 


d 2 $ 


(4.9) 


(4.10) 


(4.11) 


dagdai 

relates the internal force-like variables to the internal displacement-like variables and is 
called the internal compliance operator. This operator is completely defined once 4> is 
chosen. It is interesting to note that this operator provides information about the curvature 
of the Gibb’s potential as well as the relaxation trajectories in the associated stress space 
(Arnold [1987]). The evolution law for the internal variable is 

(4 ' 12) 


The Gibb’s and dissipation potentials are directly linked through the internal state variables. 
Clearly, this framework provides a structure in which the flow and evolutionary laws are 
fully associative and hence easily integrable. 


4.3 Generalized Viscoplastic model with Potential Structure 

Using the above framework, Arnold et al [1996] developed a multiaxial, potential based, 
fully associative, isothermal, unified viscoplastic model. This model possesses one tensorial 
internal state variable, the internal stress (ajj) that is associated with dislocation motion. 
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The forms of the Gibb’s potential (4’) and the dissipation potential (fi) are chosen as follows 
(Saleeb and Wilt [1993]): 


$ 


1 

2 


(1 + W 


G 0+1 


(4.13) 


= K 2 


1 (F) n+l 
2 fi (n + 1) 


R 

H {m + (3 + 1) 


(4.14) 


where k, //, n(> 1), ?n(> (3 + 1/2), /?, i?, and H are positive material parameters. The first 
three are associated with the flow law and the rest with the evolution equation. The scalar 
function F depends on the effective stress and G depends on the internal stress invariants. 
The McCauley brackets are such that 


(F) = F if F> 0 
= 0 otherwise 


{( G ) = G if G> Go 
= G 0 otherwise . 

This ensures use of different forms of F in different regions of the state space to account for 
stress-reversals, cyclic- loading and dynamic recovery effects. Also, G 0 is the cut-off value of 
G, which is a small constant that helps fit the experimental data and prevents singularity 
in the numerical solution when G tends to zero. 


In the spirit of von Mises and owing to the deviatoric nature of the inelastic deformation, 
only the quadratic invariants are considered in the definitions of the potentials for this 
model. Hence the definitions of F and G are as follows: 


F = 



(4-15) 


0 = 4 


where J 2 is the effective deviatoric stress invariant defined as 


J2 — 2 


(4.16) 


(4.17) 
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and J 2 is the internal deviatoric stress invariant given by 


4 


2 a ij(Hj 


and the effective deviatoric stress is 


^ ij — Sij Q>ij. 


(4.18) 


(4.19) 


The function F 1 acts like a threshold surface because no inelastic strain will occur if F < 0. 
The size of this surface is dictated by the shear strength («) and its location is determined 
by aij. Following the above framework and using these potentials (equations 4.13 and 4.14), 
the flow law is derived as 


e 1 - — £ 

^ " 2 11 13 

and the evolution equations for the internal stress are 

Q-ij = iylijki ~ (5(1 + 2j3) aijakl ') ~~ RG m a k i^j 

where Iij k i is the fourth order identity tensor. 


(4.20) 


(4.21) 


The above equation consists of two terms that represent competing mechanisms in the 
material; (1) a hardening term (associated with i/, (3) that accounts for the strengthening 
mechanisms, and (2) a recovery term (associated with R , m ) that accounts for the softening 
mechanisms and is called strain- induced recovery or dynamic recovery. These competing 
terms in the evolution equations are consistent with the assumption about the nature of 
internal mechanisms in the material (Miller [1987]; Freed et al [1991]). 


4.4 Proposed Model 

The model described in the previous section is J 2 -based, which is suitable for many metals 
like titanium (Arnold et al [1996]). It cannot, however, account for the SD phenomenon 
which is observed in aged Inconel 718. To do so, we need to generalize the above model. 
The starting point here is the choice of Gibb’s and the dissipation potentials which have the 
same form as equations 4.13 and 4.14 respectively. The threshold function is generalized 
and then the flow and evolutionary laws are derived. 


N AS A/CR— 200 1 -2 1 07 1 5 


61 



Table 4.1: Material parameters for various threshold functions 


Source 

Threshold (Yield) Function 

m 

a 

b 

c 

Mises 

V^2-k 

1/2 

0 

1/k 

0 

Drucker-Prager 

al l+v^-k 

1/2 

a/k 

1/k 

0 

Drucker 

t3 9 t2 l.6 
4 

3 

0 

1 /k 6 

-9/(4fc 6 ) 


4.4.1 Modified Threshold Function 

As mentioned before, an elastically isotropic material that exhibits a SD must be represented 
by a threshold function that depends on an odd power of I\ or J3 or a combination of 
both. Pressure sensitivity requires dependence on h while pressure insensitivity requires 
dependence on J 3 (and not /1). Here, the threshold function is generalized to include all 
the effective stress invariants as follows: 

2m | 

F = a J? mi + b j™ 1 + c J 3 3 - 1 (4.22) 

where a, 6, c, and m\ are constants to be fit to the experimental data. The invariants I\ 
and J3 are defined as 

h = &U uijt 

J3 = gSijEjfcSfcj. 

This form is convenient to use because it is an additive combination of all three effective 
stress invariants. Initially, when the internal variables are zero, the invariants depend on 
the external deviatoric stresses only, and the threshold function F is the yield function. 

There are several advantages of choosing such a form of threshold function. 

1. It is most general in that it incorporates all the three effective stress invariants. 

2. The polynomial form is convenient for differentiation which is required to derive the 
flow and evolution laws. 

3. Equation 4.22 can be easily reduced to some of the classical yield functions by choosing 
suitable values for a, b, c, and m\ as shown in Table 4.1. 
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4. For pure shear loading (7 1 =J 3 =0), the above functional form reduces to 

F — bj 2 - 1 

which has the same form as the baseline GVIPS model. This feature will be useful in 
characterizing the model, which is detailed in subsequent sections. 

Further, for better flexibility of data fitting, the above threshold function is modified as 

2m] j 

F = [a ff m +b j™' +c J 3 3 - 1. (4.23) 

Notice that threshold function in the equation 4.23 also reduces to the baseline GVIPS 
model for pure shear loading (/ 1 =J 3 =0). 

4.4.2 Flow and Evolution laws 

Function (7, which depends on the internal variables is chosen to have the same form as F 
and can be written as 

G = [al[ 2m + bJ' 2 mi + cJ^r (4-24) 

where a, 6, c, and m\ are same as those in equation 4.23. The invariants I[ and J 3 are 
defined as 

A 
4 

Using these functions in the Gibb’s and the dissipation potentials (equations 4.13 and 4.14), 
and by appropriate differentiation of the potentials, we can derive the flow and evolution 
laws. In this section, only the final formulae are given. Detailed derivation is given in 
Appendix A. 

Flow Law 

The associated flow law (equation 4.8) is used to derive the flow law. Using fi (equation 
4.14) the inelastic strain rate tensor is obtained as 

Hj = ^F"(F + Tij (4.25) 


— ota 
1 

— r> ai i a j kaki ' 
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where 


T tj = (2 a! 2 ™'- 1 - -cJ 2 J$ 


4 •- 


^’iq'Ejq- 


Equation 4.25 is the flow law for the proposed model. 


Evolution Law 


The evolution equations for the internal variables can be obtained in a similar manner using 
equation 4.12 as follows: 


Ajj 


on 

docij 


efj - ^< 5 ( 1 + 0 +*™)©. . 


u 


(4.26) 


where 


@ij = (2 al[ 


! 2m ] — 1 


- 


■/ 4 mi - 1 


)S, 


ij 


i l w m-i — 1 , ^ t/ o mi 1 

T 0^2 fly T ^^3 Q’iqO.jq 


where Oy is the internal deviatoric stress tensor. Equation 4.26 is the evolution law for the 
proposed model. 


Internal Constitutive Equation 


The internal compliance tensor Qijki is obtained by differentiating $ (equation 4.13). 

a 2 * _ k 2 81 


tyijkl r\ 

dot%j daw 

The internal constitutive rate equation is 


H dot. 


(4.27) 


*? 


oiij — LijkiAfci 


(4.28) 


where A - are displacement like variables that are conjugate to internal stress tensor (« tJ ) 
and Lijki is the internal stiffness tensor (= QjAi). Using G in equation 4.24 and differentiating 


N AS A/CR— 200 1-210715 


64 



we get 


Qijkl = (4.29) 

^ G (l+/J-m,) [(2o(2mi _ 1)7 /2m,-2 _ & + _ 1) j' 2 j' §™’ 

+6(mi - 1 )j£ m,_2 ayo« + bJp'^SikSfl 

2 

4/ + (gWl - 1) 

-y J 3 5m,_1 ((^ m i ~ VjjaicpaipSij + a kl 6ij) + (— e « e «l- 

This compliance tensor, when inverted and substituted in equation 4.28, will completely 
define the evolution of internal stress. It is seen from equation 4.27 that the internal 
compliance is not constant. It depends on the deviatoric and mean internal stresses. Thus 
equation 4.28 is highly non-linear and rate dependent. A key difference between this model 
and the baseline GVIPS model is the internal constitutive rate equation. Here, the total 
internal stress rates are related to the internal strain rate (Ajj), while in the baseline model 
the deviatoric variables are employed. 


+?£ j'l m i ' 
+ 3 3 




jf {O'iqQ'jqQ'kpQ'lp ^ a iq a jq^kl)) 


4.5 Characterization of Material Parameters 

The model now has eleven independent material parameters that need to be determined for 
aged Inconel 718 at 650°C. These parameters are summarized below. 

1. Flow law parameters (k, /i, n): k is the threshold stress in pure shear, /i is associated 
with the viscosity of the material, and n is an exponent. 

2. Hardening parameters (H, f3 ) : These are responsible for work hardening of the material 
during plastic deformation. 

3. Recovery parameters (7Z, m): These parameters are associated with the recovery mech- 
anism that competes with the hardening mechanism in the material. 

4. Threshold function parameters (a, b , c, mi): Parameter a scales the mean stress 
invariant /i while c scales J 3 , and mi is an exponent. 
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The parameter b is set to 1 /k( 2 "u) so that for a pure shear loading the threshold function, 
and hence the flow and evolution equations will reduce to the baseline GVIPS model. 

The methodology for characterizing the remaining ten material parameters is as follows: 

1. Conduct pure shear experiments on aged Inconel 718. Fit the parameters associated 
with flow and evolution laws (GVIPS parameters), k, p, n, m, /?, R and H to these 
data using optimization. Since /i=J3=0 the parameters a, c, and mi will not affect 
the prediction of GVIPS parameters. 

2. Conduct axial tension and axial compression experiments on aged Inconel 718. Use 
the GVIPS parameters (obtained from step 1) to fit the parameters a, c, and mi to 
these data. 

3. Ensure that o, b, c, and mi are such that the resulting initial threshold surface is 
convex, thereby satisfying Drucker’s stability postulate at a material point. 

These steps are necessary, in part, due to limitations with the current characterization tooLs, 
which are explained in detail in the following subsections. 

4.5.1 Parameters associated with flow and evolution 

GVIPS parameters are fit to the pure shear experimental data on Inconel 718. This is done 
using the Constitutive Material Parameter Estimator (COMPARE), a software package de- 
veloped by Saleeb et al [1998]. COMPARE was run on a NASA-GRC computer system. 
Material parameters are determined by minimizing the errors between the experimental data 
and the predicted response. It uses both optimization and constitutive model (J 2 based) 
analysis, casts the estimation as a minimum error weighted, multi-objective optimization 
problem, and then solves the optimization using a sequential quadratic programming tech- 
nique. It has three main parts that are summarized below. 
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Primal Analysis 


This module is responsible for numerically simulating the load paths used in the experi- 
ments. It is similar to nonlinear finite element code built upon a single four node plane 
stress element with material nonlinearity capability. A fully implicit backward Euler scheme 
with corresponding algorithmic (consistent) tangent stiffness matrix is used because of its 
robustness and superior stability and convergence properties. 

Sensitivity Analysis 

Sensitivity analysis involves the calculation of parameter sensitivities to the predicted re- 
sponse. The sensitivity analysis is of the direct type performed on the basis of an explicit 
recursive form associated with the above integrator. The sensitivities are naturally derived 
from exact expressions in conformity with the underlying integration scheme. The advan- 
tage of this approach is the improved computational efficiency, while the disadvantage is the 
necessary analytical derivation of sensitivities that at times become complicated. However, 
the former outweighs the accompanying complexity. 

Optimization 

The multi-objective optimization problem is formulated in COMPARE and solved using 
a sequential quadratic nonlinear programming technique (Schittkowski [1981]). Its salient 
features are: 

1. Design variable formulation that includes component synthesis, i.e. active/ passive 
design variables. During optimization, the active design variables are modified and 
the passive ones are not. 

2. General scaling of objective functions as well as design variables for numerical sim- 
plicity. 

3. Formulation of a single design optimization problem through a weighted objective 
function. 
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The weighted objective function IF j § ) for p tests is minimized to find n variables §. Here 
J r (§) is expressed as 


^(§) = E 

t=i 


mm 

A 0 


where 



(4.30) 


A 0 = '£ W M*o(§)\- 

i— 1 

W{ is the weighted parameter for the i th test such that 

= i- 

1=1 

/i(§) is the objective function for the i th test with initial value /* 0 (§). n p is the number 
of measurement stations along a load history. and Rp k are the k th components of 
response from analysis and test, respectively, n# is the number of measured components at 
a particular j th station. 


The input requirements and procedure followed for using COMPARE to estimate these 
material parameters are detailed in Chapter 6. 


4,5.2 Determination of threshold function parameters 

After the seven parameters associated with flow and evolution law are determined and 
having set b = l/« 2mi , it remains to determine the parameters a, c, and m \ . Material 
response to a general loading path that results in non-zero stress invariants I\ and J3 
is required to determine these parameters. Data from uniaxial tension and compression 
experiments (which result in non-zero I\ and J3) axe to be used for this. 

The J2 model (a = c —0) will predict identical responses in tension and compression and 
thus will not capture the SD effect. By suitably choosing the values for a, c, and rri\ it 
is possible to match the predicted responses in tension and compression with experimental 
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data. This is done iteratively. Multiple sets of parameters (a, c, mi) can be found that 
may fit the data equally well, resulting in fy J 2 , J 2 J 3 ? and I 1 J 2 J 3 models. The procedure for 
determining the proper combination of (a, c, mi) for Inconel 718 is explained later. This 
will determine whether the SD phenomenon is caused by the effect of hydrostatic pressure, 
or by influence of J 3 , or both. 

4.5.3 Convexity Requirement 

Convexity of the yield surface is a consequence of Drucker’s stability postulate. This means 
that the values of a, b , c, and mi are not completely arbitrary. Thus, for a given mi, 
though multiple combinations of a and c may fit the axial experimental data in tension and 
compression, only those combinations that satisfy the convexity requirement are acceptable. 

Outlined here is a procedure to verify convexity of any given function in three-dimensional 
space (Iyer and Lissenden [2000b]). Using this procedure, we can determine the limits that 
the convexity requirement places on the threshold function parameters. 

The basic requirement for convexity of a three-dimensional surface is that its curvature 
must be non-negative everywhere on the surface. That is, the curvature of the function 
in two mutually perpendicular directions in the tangent plane at any given point on the 
surface, must be non-negative. Mathematically, this requires that the curvature tensor be 
positive definite. 

An arbitrary vector m and the tangent plane to the threshold surface are shown in principal 
stress space in Figure 4.1 . The components of this vector in the normal and tangential 
directions to the threshold surface are also shown. Since u* is arbitrary, all vectors that lie 
in the tangent plane are given by ( 8^ — niTij)uj , where n* denotes the unit outward normal 
to F. The unit vectors in the tangent plane are Sij - riirij. Orthonormalizing this vector 
gives a set of three basis vectors, two of which lie on the tangent plane, the third being the 
unit normal vector itself. 
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Figure 4.1: Basis vector determination 


The curvature in the tangential plane is given by 
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where Hij is the Hessian of F given by 


Hij = 


d 2 F 
0(7 ‘i 0(7 j 


(4.31) 


(4.32) 


where < 7 j is the principal stress vector. U consists of two basis vectors that lie in the plane 
tangent to F in the principal stress space. 


Using the above procedure, material parameters associated with flow and evolution laws 
(«, /i, n, (3, m, f?, H) and prescribed combinations of threshold surface parameters (a, 
c, mi) will be determined. To find the right combination of (a, c, mi) that describes 
the material behavior accurately, we need to consider material response for load paths 
other than pure shear and uniaxial loadings. These are validation experiments that will 
be detailed in subsequent chapters. Uniaxial testing under hydrostatic pressure will be an 
added advantage. These will help us to confirm the importance of I\ (value of a) for this 
material. 
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Chapter 5 


Experimental Program 


This chapter starts with a brief review of relevant experimental work on inelastic defor- 
mation of metals. Next, the details of the test system and specimens used in the present 
investigation are presented. A detailed section on the design of experiments follows, which 
establishes a method to isolate the effect of each stress invariant on deformation in three 
dimensional principal stress space. The method is then simplified for axial-torsional space. 
Finally, the test matrix for the characterization tests of Inconel 718 is presented. 

Experimental evaluation is an important step in the development of multiaxial viscoplastic- 
ity models. In general, three types of experiments are necessary to support the development 
of any potential based models (Robinson [1985]) for high temperature structural materials. 
These are as follows: 

1. Exploratory tests : These tests guide the development of the theoretical framework 
and help examine the mathematical aspects of the framework. For a potential based 
framework, they help in developing the functional forms of the Gibb’s and the dissi- 
pation potentials. 

2. Characterization tests: These tests provide the required database for determining the 
specific functional forms and the parameters that represent a specific material for a 
specific temperature range. Usually these are simple tests like tension, compression 
and shear on relatively simple specimen geometries. 
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3. Validation tests : These tests provide the ultimate test of the constitutive model by 
comparing actual structural component response with analytical predictions based on 
the proposed model. For most cases, simple specimen geometries are used and the 
loadings are complex enough to rigorously test the model. 

Tensile and compressive tests with superimposed hydrostatic pressure provide information 
about the significance of I\ in describing the flow and evolution laws for the material. In 
the present investigation, these tests axe intended to be used qualitatively to validate the 
models that include I\ in the threshold function. 


5.1 Review of experiments 

Plasticity experiments have been conducted on metals to study their yield and flow behavior. 
Volume expansion during inelastic deformation was the emphasis in many investigations. 
Both uniaxial and multiaxial loadings have been reported. Such experiments date back to 
Lode [1926] who tested thin walled metal tubes in tension and internal pressure. There 
was a discrepancy between experimental data and the Mises yield criterion, which led 
Taylor and Quinney [1931] to further investigate the effect of biaxial (axial and torsional) 
loadings on tubular specimens of copper, aluminum and steel. In general, their work showed 
good agreement with Mises yield criterion for these materials. They also performed large 
deformation uniaxial tensile tests on these materials and monitored the changes in internal 
volume of the specimens. They found large changes in internal volume that could not be 
accounted for by small density changes and attributed it to the anisotropy induced during 
inelastic deformation. 

Most biaxial testing was done using proportional loading (Phillips et al [1957,1961], Findley 
et al [1962], Michno and Findley [1976]). Biaxial experiments on metals at high tempera- 
tures (> 500° C) were difficult to perform owing to the complexity in experimentation and 
the associated difficulty in accurate measurement of strains. 

Dependence of inelastic deformation on hydrostatic pressure has been another topic of 
investigation for metals. Contrary to the experimental findings of Bridgman [1952], some 
metals have exhibited dependence of flow on hydrostatic pressure. Drucker [1973] pointed 
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out that the SD phenomenon can result as a consequence of pressure dependent yielding. 
Spitzig et al [1975] have shown that several high strength steels (tempered 4310 and 4330) 
that exhibited SD effect indeed showed pressure dependent yielding. The volume expansion 
observed in these materials was in accord with the expected increase in dislocation density 
during deformation, but did not agree with the normality rule predictions. 

Another type of plasticity experiment is the determination of initial and subsequent yield 
surfaces. Comprehensive reviews of these experiments have been provided by Michno and 
Findley [1976] and Hecker [1976]. Most of these experiments involved axial torsional loading 
on thin walled tubular specimens. The pros and cons of using strain and stress controlled 
loadings were investigated by some researchers. Phillips and Lu [1984] used stress and 
strain controlled loadings to determine yield surfaces for pure aluminum and found no ap- 
preciable difference between the two. On the other hand, some researchers preferred strain 
controlled loading to stress controlled loading and argued that the former leads to more 
accurate results(Wu and Yeh [1971], Ellis et al [1983]). Dependence of yielding on tex- 
ture was determined by Althoff and Wincierz [1972] by experiments on textured brass and 
aluminum. Subsequent yield surfaces were also a focal point in many experimental investiga- 
tions (Nagdhi et al [1958], Phillips and Tang [1972], Williams and Svensson [1970,1971]). In 
general, the inelastic strain rate vector was found to be normal to the yield surface (Michno 
and Findley [1974], Phillips and Moon [1977], Khan and Wang [1993]). This suggests that 
an associated flow law is a good approximation for most metals. 

The shape of flow surfaces is strongly dependent on the flow definition (chapter 3). There 
has been some effort in the determination of flow surfaces, both surfaces of constant inelastic 
power (SCIPs) (Clinard and Lacombe [1988]) and surfaces of constant inelastic strain rate 
(SCISRs) (Battiste and Ball [1986]). These flow surfaces are either determined directly 
during experiments or by using post experiment data reduction techniques. For direct 
determination, the inelastic strain rate must be calculated in real time (Ellis and Robinson 
[1985], Battiste and Ball [1986], Lissenden et al [1997a]). 

Gil [1999b] established an experimental program to determine yield and flow surfaces for 
solutioned and aged Inconel 718. Using axial-torsional loading SCISRs and SCIPs in axial- 
shear stress plane were constructed (for temperatures between 23°C and 649°C). He found 
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that the initial yield surfaces for solutioned Inconel 718 at low temperatures fit the Mises 
criterion. However, the yield loci for solutioned Inconel 718 at high temperature (> 350°C) 
and the yield loci for aged alloy for all the temperatures investigated, showed an eccentricity 
toward compressive stress direction (see Figure 1.1(b)). This is the SD phenomenon in 
Inconel 718 and is the topic of present investigation. 

In general, plasticity experiments on metals have dealt with determining yield surfaces (ini- 
tial and subsequent) and finding the effect of loading (type of control, rate) on hardening 
behavior. The effect of hydrostatic pressure on yield and flow is another topic of interest. 
The present work involves loadings, both uniaxial and biaxial (axial- torsional), deep in the 
inelastic region for aged Inconel 718 at 650° C. Elevated temperature experiments involv- 
ing multiaxial loading are difficult to perform for two reasons; (1) difficulty in accurate 
measurement of strains at elevated temperature and (2) the need for such experiments was 
not felt because no systematic effort was made to determine the individual effects of stress 
invariants on the inelastic flow of metals. Of course, in the present investigation there is 
such a requirement in order to capture the SD phenomenon. A comprehensive technique to 
capture the SD effect by investigating the effect of all the three stress invariants on inelastic 
deformation is not reported to date. 


5.2 Test equipment and Specimen details 

An MTS axial-torsional test system was used for the uniaxial and biaxial experiments. 
A gas based high pressure deformation apparatus was used for testing under hydrostatic 
pressure. The specimens used for testing have different geometries for each of these systems. 

5.2.1 MTS Axial- Torsional test system 

The MTS test system (at NASA-Glenn Research Center) is a servo-hydraulic test machine 
having an axial load capacity of 222500 N and a torque capacity of 2260 N-m. The specimens 
are gripped by hydraulically actuated grips. Shown in Figure 5.1 is the complete test system 
with all the required accessories. The top grip of the load frame is attached to an axial- 
torsional load cell that in turn is attached to the cross head, which remains fixed during 
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Figure 5.1: MTS Axial-Torsional test system at NASA Glenn Research Center 


a test. The bottom grip is attached to an actuator capable of independent rotational and 
vertical translational motions. The rotation of the actuator can be controlled in a closed 
loop system by the angle of rotation, torque and shear strain while the vertical motion is 
controlled by either displacement, load or axial strain. Kalluri and Bonacuse [1990] provide 
additional details regarding the biaxial test machine. 

The test machine is equipped with an adjustable coil (figure 5.2), 50 kW audio frequency, 
induction heating system capable of generating specimen temperatures up to 800°C. For 
this reason, the specimen grips axe water cooled. The temperature control is done by 
one thermocouple spot welded to the the specimen at the gage section. The temperature 
distribution in the gage section is determined by thermocouples which are spot welded to 
the outer surface of the specimen in the gage section. The temperature variation is limited 
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Figure 5.2: Close up showing the specimen, heating coil and extensometer 


to ±1% of the target temperature. An enclosure around the test machine limits the effects 
of air currents. Axial and torsional stress and strain data are saved electronically. 


Strain measurement 

At high temperatures like 650° C it is not practical to use strain gages. On the other hand, 
strain measurement with good resolution is required to accurately determine the initiation 
of yield. At high temperatures extensometers are preferred to strain gages. The common 
factors that make accurate strain measurement (using extensometer) at high temperatures 
difficult are (1) presence of electronic noise and (2) coupling between axial and torsional 
strain components. 

A multiaxial extensometer (figure 5.2) that is capable of axial and torsional strain measure- 
ments over a wide range of temperatures is used. The extensometer contains two alumina 
rods spaced 25 mm apart. These are in contact with the specimen by means of indentations 
on the specimen and spring loading provided by the mounting fixture. The top rod is free 
to move only in the axial direction while the bottom rod is free to move only in the cir- 
cumferential direction. Axial displacement S and angle of twist 0 are output voltages. The 
axial strain is en = S/l 0 and the shear strain is 712 = r o 0/l o , where r 0 is the outer radius 
of the specimen and l Q is the gage length (25 mm) of the extensometer. 
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servohvdraitlic 

actuator 


Figure 5.3: Gas based high pressure deformation apparatus 

Test machine control 

The MTS system is controlled by a microcomputer equipped with a digital-analog (D/A) 
converter that provides independent controls over the axial and rotational motions of the 
actuator. A 16-bit analog-digital (A/D) converter is used to acquire load, torque and 
extensometer (axial and torsional) data. The D/A and A/D hardware is commanded 100 
times per second by software that is customized using a FORTRAN program. Two tests 
were performed using stress control. All other tests in this investigation were performed in 
strain control. 

5.2.2 Pressure test equipment 

A schematic diagram for the gas based high pressure deformation apparatus (at Case West- 
ern Reserve University (CWRU)) is shown in Figure 5.3. It utilizes a pressure intensifier 
to generate pressure that is contained within the multi-walled pressure vessel. The volume 
of the pressurized gas is kept as low as possible because of the danger associated with the 
stored energy. Pressure is monitored using manganin coil pressure gage that is exposed to 
the high pressure environment. These coils are used because of the highly reproducible and 
linear manner in which the coil resistance changes with applied pressure. 
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Figure 5.4: Dimensional details of tubular specimens (all dimensions in mm) 

The specimen is first inserted in the load train assembly present in the pressure vessel 
before gas pressurization and then tension or compression is applied at a desired level of 
superimposed hydrostatic pressure. It is important to continuously monitor the pressure 
during testing. It is required to accurately monitor the load and displacement during 
deformation under pressure. For more details on the test system and test procedures, the 
reader is referred to Lewandowski et al [1998]. 

5.2.3 Specimen details 

Specimens used in the MTS system were thin-walled tubes, designed to achieve a plane state 
of stress in the gage section. Also, the dimensions are such that loading the gripped ends of 
the specimen results in a uniform stress distribution in the gage section. The dimensional 
details of the specimen are given in Figure 5.4. 


5.3 Design of Experiments 

We have introduced the importance of stress invariants other than J 2 to represent yield and 
flow behavior for a class of materials that exhibit a SD effect. These invariants are brought 
into the formulation by simply including them in the threshold function and then deriving 
the flow and evolutionary laws based on it. Although the chosen form of the threshold 
function in equation 4.22 is justified for its simplicity and easy reducibility to well known 
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forms, it is only one of the many possibilities. The appropriate choice of the form may 
depend on the material under investigation. 

A traditional approach for metals would be to conduct hydrostatic tests as suggested by 
Bridgman [ 1952 ]. For pure hydrostatic state of stress (cr), I\ = 3 a and J2 ~ J3 = 0 - Hence, 
such a test will help determine whether I{ affects flow behavior and causes permanent 
volume change during inelastic deformation. If no permanent deformation is observed under 
hydrostatic stress, the threshold function can be based on deviatoric invariants alone. Unlike 
the threshold function for a J2 material, which can be quantified by pure tensile loading, 
a simple uniaxial test will be insufficient to quantify a general threshold function. This 
is because the general yield function has three unknown coefficients, the determination of 
which requires at least three experiments. 

5.3.1 Isolating the effect of invariants 

For a general case where I\ and/or J3 must be considered, a multiaxial test program needs 
to be adopted to quantify their relative importance in the definition of yield and flow. This 
can be accomplished by following stress trajectories that have only one stress invariant that 
is changing, or better yet, only one nonzero stress invariant. Some simple loadings (external) 
such as hydrostatic pressure (/1 7^ 0, J2 = J3 = 0) and pure shear (J2 7^ 0, I\ = J3 = 0) 
can help describe the initiation of inelasticity. The internal state of the material does not 
change during elastic deformation and hence I\ = /1, J2 = J 2 , and J3 = Js- During inelastic 
deformation, however, the external loading that causes pure shear or hydrostatic tension are 
not immediately known. This is because the invariants also depend on the internal stresses 
which are nonzero during inelastic deformation giving I\ / Ji, J2 7^ J2, and J3 7^ ^3- 

For a general loading condition, stress paths having only one changing effective stress in- 
variant can be obtained by following guiding vectors in the effective stress plane determined 
by the cross product of the gradient of the other two invariants. Mathematically, this can 
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be expressed as 


VJ2 x VJ3 : /1 changing 

VJ3 x V/i : J2 changing ( 5 . 1 ) 

V/i x VJ2 : J3 changing 

where Ji, J2 and J3 are effective stress invariants defined previously. 

Stress paths obtained by the above equations are greatly influenced by the starting point, 
which is not required to be the zero stress condition. Effective invariants also depend on 
the internal stresses, which are neither measurable nor controllable variables. Therefore, 
computer modeling is necessary to determine what loading to apply in order to follow 
the guiding vectors in equation 5 . 1 . An additional complexity arises for non-proportional 
loading for which the applied stress and effective stress directions are generally not the 
same. 

5.3.2 Stress trajectories in three dimensional stress space 

Using equations 5 . 1 , a number of stress trajectories can be derived depending on the starting 
point in the stress space. Figure 5 . 5 (a) shows two stress paths in the principal effective 
stress space that have both J2 and J3 constant, while Figure 5 . 5 (b) shows the change in I\ 
along the same stress paths. Note that path AB is hydrostatic tension while path CD is 
not. However, both paths are straight and have the same direction. Similarly, the case for 
varying J2 with constant I\ and J3 is shown in Figure 5 . 6 . Again, the importance of the 
starting point is clearly seen, path AB starts very near the origin and is pure shear while 
path CD is not. Both AB and CD are straight and parallel but the increase of J2 along CD 
is larger than that along AB. 

Finally, paths for varying J3 with constant I\ and J2 are shown in Figure 5 . 7 . They are 
neither straight nor do they intersect the origin. Both paths have a local extreme which is 
expected to cause a distinct change in the inelastic response as the change from loading to 
unloading occurs. 
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Figure 5.5: Effective stress paths with I\ changing - J 2 and J 3 constant 



N AS A/CR— 200 1-210715 


81 








(MPa ) 


1 2 (Mr a) .500 .500 - 400 *1 < Mra ' 



Dist. along path (MPa) 


Point 

A 

r l 

-500 

r 2 

400 

r 3 

250 

B 

-94 

-344 

588 

C 

-245 

133 

512 

D 

453 

-285 

231 

AB 

b = 

J 2 - 

150 MPa 

2.325 x 10 5 (MPa/ 


L = 400 MPa 

CD ~ k 2 

J 2 = 1.432 x 10 5 (MPa) 


Figure 5.7: Effective stress paths with J 3 changing - I\ and J 2 constant 
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Experimentation using three-dimensional stress states require complex laboratory equip- 
ment such as that developed by Calloch and Marquis [1999]. Hence, a simpler alternative 
approach would be to restrict the stress trajectory to a plane stress. Plane stress loading 
can be readily applied to a thin-walled tube by axial force, torque, and internal pressure. 
For a plane state of stress however, the invariants are related by 

h = (5.2) 

resulting in only two of the invariants being independent. Thus a different approach than 
that for three-dimensional space is needed. 

5,3. 3 Plane stress experiments 

Since J 2 has the primary influence on the inelastic flow, it can be treated differently relative 
to the other two invariants. It is thus proposed that the change in J 2 be the same in each 
of the two tests, one test having a varying I\ and a constant J 3 and the other test having a 
varying J 3 and a constant I \ . These types of tests appear to be possible using axial-torsional 
loading of a thin walled tubular specimen. 

Simple load paths in the axial-torsional stress space that satisfy the above condition are 
shown in Figure 5.8 and Figure 5.9. If a pure tensile stress state is chosen as a starting 
point for shear loading as in Figure 5.8, the result is a constant I\ and a changing J 3 . If a 
combined compression-shear stress state is chosen as a starting point and the compressive 
stress is reduced while the shear stress is increased as in Figure 5.9, J 3 remains constant 
and I\ increases. Though I\ approaches zero, this load path corresponds to loading since 
the J 2 rate is positive in the sense of increasing inelasticity. 

The intent of these tests is to facilitate comparison with pure shear tests conducted at an 
identical J 2 rate. Differences in the inelastic strain response could then be attributed to 
J 3 for loading in Figure 5.8 and to I\ in Figure 5.9. This can be used to determine an 
appropriate weighting of the material parameters a and c in equation 4.22. 

Since proportional loading is easier to apply than non-proportional loading, it is worthwhile 
to consider the proportional load paths shown in Figure 5.10. The loading rate can be 
adjusted so that the J 2 rate is the same for any loading direction, Ai. For values of Ai of 
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Figure 5.9: Stress path with J 3 constant, J 2 rate constant and I\ changing 
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Figure 5.10: Effective stress invariants for two proportional load paths 


0.0 and 0.2816, the changes in J 3 happen to be almost identical while the change in Ii is 
not. Thus any difference in the inelastic response could be attributed to I\. 

The objective of the various experimental possibilities described above is to determine what 
effect each of the three stress invariants have on inelastic deformation. These paths can 
only be determined by using the viscoplastic model to predict inelastic strain and searching 
for paths that extremize the difference between these forms. It is probable that these paths 
will require non-proportional loading and creep or relaxation periods. 
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5.4 Proposed experiments for the current model 


First, characterization is done using pure shear tests that are similar to load path AB in 
Figure 5.6(b). Next, axial test results are used to determine the threshold surface parameters 
a and c. The intention here is to follow load paths that result in varying 7i, J 2 and J 3 , and 
hence are not similar to any of the paths in Figures 5.5 to 5.10, where the idea was to keep 
two invariants constant and vary the third. 

Validation experiments are biaxial tests using, both proportional and non-proportional load- 
ing. The non-proportional loading is similar to the path in Figure 5.8 while the proportional 
loading is similar to the load paths in Figure 5.10. Tension and compression tests with super- 
imposed hydrostatic pressure complement these validation experiments since they directly 
show the effect of 7i on material inelasticity. Test matrices are developed for characteriza- 
tion and validation of the proposed viscoplastic model in the next section. 


5.5 Test Matrices 

The most important and often times the most difficult aspect of modeling at elevated 
temperature is obtaining the required material parameters. The associated difficulty stems 
not only from the variety in the mathematical forms of the threshold function, but also from 
the fact that multiple sets of material parameters can correlate experimental data equally 
well, for a given load path. In order to arrive at the proper set of material parameters, it is 
therefore crucial to choose an appropriate set of experiments that bracket a wide range of 
values for the variables under consideration (e.g. loading rate, load levels, control mode). 

First, the test matrix for the pure shear characterization is developed. This is followed by 
tests in tension and compression. Finally, test matrix for axial testing under hydrostatic 
pressure is also given. All tests axe conducted at 650°C except for the pressure tests, which 
are at room temperature. 

Pure shear characterization 

Since we are trying to determine seven material parameters, different types of experiments 
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Table 5.. 

L: Pure shear characterization experiments 

Name 

Control 

Loading rate 

ts(sec) 

tf(sec) 

End level 

IN7 

Strain 

1745 microstrain/s 

10.8 

43200 

0.0188 

IN20 

Strain 

17.4 microstrain/s 

1500 

44700 

0.0261 

IN21 

Stress 

25.1 MPa/s 

22.9 

1282 

552 MPa 

IN22 

Stress 

25.1 MPa/s 

19.5 

43200 

482 MPa 



Figure 5.11: Schematic for shear characterization tests 


are typically required for proper estimation of these parameters. On this basis, four tests 
(summarized in Table 5.1) were conducted. 

The strain controlled tests were at different loading rates (two orders of magnitude) and 
have the same end strain level. On the other hand, the stress controlled tests were at the 
same loading rate but have different end stress levels. All the hold times except for IN21 
were approximately 12 hours. A schematic diagram for these tests is shown in Figure 5.11. 

Axial tests 

Axial tests that will be used to characterize the material parameters a, c and mi are 
summarized in Table 5.2. 

Testing under hydrostatic pressure 

The intention of these tests is primarily to get a qualitative estimate of the contribution of 
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Table 5.2: Axial characterization experiments 


Type 

Control 

End level 

Time 

Tension 

Strain 

0.02 

100 s 

Compression 

Strain 

-0.02 

100 s 


/ [ in the threshold function. Hence, only a few tests are planned which are either tensile 
or compressive tests under three different values of hydrostatic pressure. These tests are 
planned to be done at room temperature. The pressures are in the range; 0.1 to 420 MPa. 

After obtaining the material parameters using shear and axial characterization, validation 
of the model is done using biaxial experiments, which are developed in the next chapter. 
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Chapter 6 


Results and Discussion 


We have now developed the formulation for a unified viscoplastic model using a generalized 
threshold function and identified experiments that are required to determine material pa- 
rameters in the model. The results of the characterization and validation tests are presented 
in this chapter. The procedure for determination of the material parameters is explained. 
Once the parameters are determined, biaxial validation tests are compared with model 
predictions. Results from the tension and compression tests under hydrostatic pressure il- 
lustrate the importance of the first stress invariant in the description of inelasticity. Finally, 
the results are discussed in detail and their implications on the model are described. 


6.1 Determination of material parameters 

The need for a unified viscoplastic model capable of capturing the SD effect has been estab- 
lished and a generalized theory has been developed. Additionally, an experimental program 
to characterize the material parameters in the model was developed. In this section, the 
experimental results are presented starting with pure shear experiments. Using these data, 
the parameters associated with GVIPS are optimized using the program, Constitutive Ma- 
terial Parameter Estimator ( COMPARE , Saleeb et al [1998]) and the correlated responses 
are presented. COMPARE optimizes the parameters associated with flow and evolution (k, 
n, /i, m, /?, P, H). Determination of threshold function parameters (a, 6, c, mi) is done using 
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axial (tensile and compressive) experimental data. Once all the parameters are determined, 
the models are put through a rigorous test to verify convexity of the threshold functions. 

6.1.1 Experimental results 

The pure shear experiments described in Table 5.1 were carried out on aged Inconel 718 
at 650° C. A wide range of loading rates and end levels are used in the characterization 
experiments and hence we can expect to obtain a good estimate of the material parameters. 

Results of the pure shear experiments are presented in Figure 6.1 to Figure 6.3. Stress-strain 
response for monotonic increasing shear strain at two different rates is shown in Figure 6.1. 
The loading rate for IN 7 (1740 microstrain/sec) is two orders of magnitude higher than 
that for IN20 (17.4 microstrain/sec) resulting in less inelastic deformation in IN7. Hence, 
the maximum stress level for IN7, in spite of its lower strain level, is higher than IN20. 
Figure 6.2 shows stress relaxation over 12 hours after loading these two specimens. Stress 
relaxation is higher for IN7 due to the smaller amount of inelastic deformation that occured 
during loading relative to IN20. Creep test results are shown in Figure 6.3. Creep strains for 
IN21 are larger than IN22 because the constant shear stress was higher; 552 MPa compared 
to 482 MPa (for the same loading rate of 25.1 MPa/sec). 

Strain controlled loading is chosen for further uniaxial and biaxial experiments. For this 
reason, only the strain controlled test results from IN7 and IN20 are used for characterizing 
the material. 

6.1.2 Optimization using COMPARE 

The size of the data set was reduced to about 12-15 points for each test in order to run 
the optimization program in a reasonable amount of time (less than approximately two 
hours). Care was taken in order not to lose essential information provided by these tests. 
The reduced test data points are shown by circles in Figures 6.4 and 6.5. 
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Figure 6.1: Shear stress-strain response 
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Shear Stress (MPa) 



Figure 6.2: Shear stress relaxation 
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Table 6.1: Experiments used for shear characterization 


Name 

Max. Strain 

Time to max. strain (sec) 

Hold time (sec) 

Weights 

IN7 

0.0188 

10.8 

- 

0.25 

IN20 

0.025985 

1499.2 

- 

0.25 

IN7 

0.0188 

10.8 

43200 

0.25 

IN20 

0.025985 

1499.2 

43200 

0.25 


Table 6.2: Elastic constants for aged In conel 718 at 650°C 


E 

165360 MPa 

V 

0.297 

G 

63732.5 MPa 


The choice and the number of points used in the test data plays an important role in 
optimizing the material parameters to fit the data. As can be clearly seen, more points are 
used in the regions of changing slopes. These regions are typically transitions from elastic 
to inelastic regimes or from monotonic loading to hold. 

COMPARE was used to fit the parameters associated with the flow law, hardening, and 
recovery parameters using these four sets of experimental data. The input data are sum- 
marized in Tables 6.1 and 6.2. 

An intermediate step was adopted to facilitate optimization of the seven material parameters 
(rc, n, p, m, /3, R, H). First, only the loading portions of IN7 and IN20 were used with equal 
weights (0.5) to evaluate all of the parameters (set I, Table 6.3). Next, the complete data 
sets from IN7 and IN20 (including relaxation) were used to determine another set (set II) 
of par am eters (Table 6.4). When all four tests are used together, the parameters associated 
with flow law («, p, n) were kept close to set I and those associated with evolution (m, 3, 
R, H) are kept close to set II. Parameters /i, R and H were kept within ± three decades and 
parameters k, n, m and (3 were kept within ±30%. The initial values and bounds in Table 
6.5 reflect this approach. 

COMPARE successively updated the objective function until it reached a value of 1.0 and 
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Table 6.3: Optimized values of material parameters (Set I) 


Parameter 

Final value 

k (MPa) 

78.0 

n 

7.4 

H (MPa-s) 

7.0 xlO 6 

m 

6.8 

(3 

3.8 

R (1/s) 

3.0xl0” 13 

H (MPa) 

4.2xl0 10 


Table 6.4: Optimized values of material parameters (Set II) 


Parameter 

Final value 

k (MPa) 

270.0 

n 

10.2 

H (MPa-s) 

8.0 xlO 8 

m 

12.0 

P 

3.05 

R (1/s) 

6.0xl0 -13 

H (MPa) 

2.0xl0 9 


Table 6.5: Material parameters to be optimized 


Parameter 

Initial value 

Lower bound 

Optimized Value 

Upper bound 

«(MPa) 

78.0 

60.0 

93.447 

100.0 

n 

7.4 

5.0 

9.445 

10.0 

fj, (MPa-s) 

7.0 xlO 6 

7.0 xlO 3 

6.544 xlO 9 

7.0 xlO 9 

m 

12.0 

8.0 

9.9 

16.0 

(3 

3.05 

2.0 

3.445 

4.0 

R (1/s) 

6.0xl0 -13 

6.0x 10~ 16 

5.787xl0 -15 

6.0x 10 _1 ° 

H (MPa) 

2.0xl0 9 

2. Ox 10 6 

9.332xl0 8 

2.0xl0 12 
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the optimized material parameters are given in Table 6.5. To complete the optimization, 
COMPARE was required to be restarted a few times (typically 3 to 4). Each time the 
program is restarted, a new (updated) set of initial parameters is used which ensures that 
the optimization does not stop at a local minimum. The resulting correlation to the experi- 
mental data is shown by solid lines in Figure 6.4 and Figure 6.5. There is a good correlation 
between the model response and the experimental data. It should be noted that predictions 
from COMPARE using the parameters in Table 6.5 for axial loading results in identical 
stresses in tension (figure 6.6) and compression except for the sign. This is because the 
threshold function in the COMPARE (GVIPS formulation) depends only on J 2 . 

Figure 6.6 shows that the flow stresses predicted by COMPARE are higher than the exper- 
imental values for tension and compression. Positive values of either a or c or both result 
in a decrease in tensile flow stresses and an increase in compressive flow stresses relative to 
the J ‘2 model (a =c =0). Hence, it is not possible to start with the optimized parameters 
in Table 6.5 and get good correlation with tensile and compressive experiments. In order 
to successfully capture the SD effect, it is required to obtain an axial prediction that is 
between the tensile and compressive experimental data. Starting with such a J 2 model ( a= 
c= 0) we can then capture the SD effect by introducing positive values for a or c or both. 
To do this, we used an alternative approach to determine the material parameters. 

One reason for the over-prediction of flow stresses in tension and compression could be a 
need for more variety in the experimental data. An attempt to use the stress controlled 
test data (IN21 and IN22) in addition to those in Table 6.1 did not result in any better 
correlation with experiments. Even differential weighting of the tests for optimization gave 
no improvement in prediction. One approach would be to add experiments with different 
loading rates, hold times, end levels, modes of control (stress or strain) and types of loading 
(single step or multi-step), in the characterization. However, we restrict additional shear 
experiments to only one (IN4) and choose an alternative way to characterize the GVIPS 
parameters. 
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Figure 6.4: COMPARE fit for Shear loading: (a)IN7 (b)IN20 
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Figure 6.6: Prediction for axial loading: Comparison with IN6 and IN2 
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Table 6.6: Revised experiments used for characterization 


Type 

Name 

Max. Engineering Strain 

Loading time (sec) 

Weights 

Shear 

IN7 

0.0188 

10.8 

0.25 

Shear 

IN4 

0.024 

100.0 

0.25 

Axial 

IN6 

0.02 

100.0 

0.4 

Axial 

IN2 

-0.02 

100.0 

0.1 


6.1.3 Alternative approach 

A good way to start is to look at the parameters, which had been fit to the experiments 
purely by a mathematical technique, from a physical viewpoint. The flow law parameters k 
(93.447 MPa) and n (9.445) from Table 6.5 seem to be physically unrealistic, k represents the 
initial shear yield strength (equation 4.15) and hence should be close to the experimentally 
determined value of 220 MPa. Parameter n represents rate sensitivity of the material. A 
large n (9.445) gives large rate-sensitivity, but this is not observed in the experiments. Also, 
the hardening parameter H (9.332xl0 8 MPa), which represents the internal stiffness of the 
material, should be less than the external stiffness (165360 MPa). 

More reasonable values of k, ti and 77, used in the optimization, resulted in better comparison 
with axial test data. However, the predicted axial stress was still not between the tension 
and compression test data. At this point, it was felt that using different experimental data 
for the characterization would give better axial response predictions. First, test IN20 was 
removed from the characterization because its loading rate (17 microstrain/s) was much 
slower than what is planned for the present investigation (200-300 microstrain/s). Also, 
since we are not looking at long term relaxation behavior at this point, only the loading 
part of IN7 was considered. Another strain controlled shear loading test (IN4 shown in 
Figure 6.7) was added to the data because its loading rate (240 microstrain/sec) is in 
the 200-300 microstrain/s range. Further, tensile and compressive test data (IN6 and IN2, 
respectively) are incorporated in the characterization with appropriate weighting. The tests 
used in the revised COMPARE optimization are summarized in Table 6.6. The weights for 
IN6 and IN2 were developed after a few iterations. 
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Shear Si 



Engineering Shear Strain 

Figure 6.7: Shear stress-strain response (IN4) 
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Table 6.7: Optimized values using re vised tests 


Parameter 

Final value 

k (MPa) 

234.7 

n 

2.0 

\x (MPa-s) 

1.7681 xlO 5 

m 

7.0 

P 

3.652 

R (1/s) 

l.OxlO -12 

H (MPa) 

1.966xl0 4 


The corresponding optimized material parameters are given in Table 6.7 and the correlated 
shear responses are shown in Figure 6.8. The correlation for IN7 is excellent and that for 
IN4 is reasonable. In addition, the presence of data for IN6 and IN2 in the characterization 
helped to obtain an axial response that is between the experimental tensile and compressive 
data (Figure 6.9). 

Such a correlation is acceptable because of the following: 

1. The value of k is reasonably close to the shear strength (220 MPa). 

2. The lower rate sensitivity observed in experiments is obtained with a lower value of 
the exponent n (2.0). 

3. The internal stiffness parameter H is less than the elastic modulus. 

Incorporation of I\ and J% through positive values of a and c will result in more accurate 
predictions in tension and compression and thus will represent the SD effect. Of course, from 
this point onwards it is required to use the formulation developed in Chapter 4, which uses 
a general yield function (we cannot use COMPARE anymore). The detailed formulation of 
the constitutive equations and simplifications for various cases are presented in Appendix A. 
A FORTRAN program was written to calculate the material response using this generalized 
formulation and is listed in Appendix B. 
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Figure 6.8: Revised predictions for shear loading: (a)IN7 (b)IN4 
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Figure 6.9: Revised correlation for axial loading: Comparison with IN6 and IN2 
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6.1-4 Axial Characterization 


The aim of the axial characterization is to determine the threshold surface parameters a, 6, 
c and mi, which along with the parameters in Table 6.7 will form the complete parameter 
set required for the model. First, a value of mi is arbitrarily chosen. The parameter b is 
equivalent to 1 jn lmx . By doing so, the general threshold function (equation 4.23) reduces 
to equation 4.15 of the baseline GVIPS model for pure shear loading (/i=J3=0). This 
correspondence is necessary because the parameters in Table 6.7 are based on the GVIPS 
model (equation 4.15) but are also used in the general formulation. Two values for mi (1.0 
and 0.6) are investigated. Values of mi higher than 1.0 do not significantly improve the 
axial predictions. Moreover, they result in very small numbers for a and c, which cause 
numerical difficulties. The values of a and c are varied iteratively until the tensile and 
compressive flow stresses predicted by the program given in Appendix B correlate well with 
the experimental data. 

For each value of mi, three combinations of invariants were investigated. These three sets 
can be classified as 

1. 1 1 J 2 model; c = 0 and a ^ 0 

2. J 2 J 3 model; a = 0 and c ^ 0 

3. I 1 J 2 J 3 model; a / 0 and c / 0. 

The idea of choosing I\ J 2 , J2J3? and I1J2J3 models is to clearly identify the influence of 
each invariant on the inelastic behavior of the material. The threshold surface parameters 
are summarized in Table 6.8. 

The correlation of the models in Table 6.8 with experimental data is shown in Figure 6.10 
and Figure 6.11. Each of these models gives a fairly good correlation of the flow stresses in 
tension and compression. The point to note here is that all of these models do an excellent 
job in the prediction of SD especially in the regions of high inelasticity (>1% strain). This 
was our primary objective. 

Before proceeding on to the validation of the models in Table 6.8, we will compare the 
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Table 6.8: Threshold surface parameters 


Name 

m i 

6(l/MPa 2mi ) 

a(l/MPa 2mi ) 

c(l/MPa 2mi ) 

h 

1.0 

1.8154x 10 -5 

0 

0 

hh 

1.0 

1.8154x 10" 5 

5.3xl0 -7 

0 

hh 

1.0 

1.8154xl0“ 5 

0 

4.5x 10~ 6 

hhh 

1.0 

1.8154xl0" 5 

2.6xl0 -7 

2.2xl0“ 6 

h 

0.6 

1.43xl0 -3 

0 

0 

hh 

0.6 

1.43 xlO -3 

4.5 xlO -5 

0 

hh 

0.6 

1.43xl0 -3 

0 

1.5xl0 -4 

hhh 

0.6 

1.43xl0“ 3 

2.2 xlO -5 

7.0 xlO -5 


initial threshold surfaces with experiment and verify their convexity. 

6.1.5 Comparison with experimental threshold surface 

The initial threshold surface for aged Inconel 718 at 650°C in the axial-shear plane was 
determined experimentally by Gil [1999b] using an offset strain definition of 30 microstrain. 
This is shown by circles in Figure 6.12. The initial threshold surface data is regressed 
(Miller et al [1990]) to a I\ function to get the experimental threshold surface shown by 
+ symbols (mi — 1.0, a = 3.838xl0“ 6 , b = 1.8263xl0” 5 ). Also plotted are the threshold 
surfaces predicted by the model using the parameters shown in Table 6.8 for mi = 1.0. 
Comparison between the predicted and experimental surfaces clearly shows that the models 
under-predict the tension-compression yield asymmetry. 

In aged Inconel 718, the SD starts with a high value, then rapidly decreases and attains 
a constant value as inelasticity occurs (Iyer and Lissenden [2000b]). We are interested in 
material behavior as it evolves during inelastic deformation. Thus, we ignore the under- 
prediction of the SD at flow- initiation as long as the prediction is good at moderate to large 
inelastic strains. With this in mind, we choose the threshold parameter values in Table 6.8 
and ignore the parameter set obtained by regression to the experimental yield surface data. 
Moreover, use of parameter set obtained by regression under-predicted the flow stresses in 
both tension and compression. 
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Axial Stress (MPa) 



Axial Strain 


Figure 6.10: Tensile and Compressive loading: Models J2, I\ J2, J2J3, h J2J3 (mi=1.0) and 
Experiments (IN6, IN2) 

6.2 Convexity of initial threshold surface 

Convexity of the initial threshold surface is a requirement that is a consequence of Drucker’s 
stability postulate. Figure 6.12 shows the initial threshold surfaces for the parameters corre- 
sponding to mi =1.0 in Table 6.8. For mi =0.6 the surfaces are identical. As the contribution 
of J3 increases (value of c relative to a) there appear to be non-convexities appearing at 
the regions of pure shear. It is difficult to check the convexity of these functions visually. 
Hence, it needs further attention and a rigorous check of convexity for each parameter set 
in Table 6.8 is provided in this section. 
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Axial Stress (MPa) 



Axial Strain 


Figure 6.11: Tensile and Compressive loading : Models J 2 , I\ J 2 , «/2«/3> I 1 J 2 JZ (mi =0.6) 
and Experiments (IN6, IN2) 








6.2.1 Mathematical Implementation 


The objective here is to verify the convexity of the threshold function in three dimensional 
stress space. We start with the function 

F = a/j mi + bJ + cJ 3 2m,/3 - 1 (6.1) 


and follow the general methodology to check convexity, outlined in chapter 4. For ease of 
computation the principal stresses (cti, ct 2 , a 3 ) are expressed in terms of the octahedral 
normal stress ( <r oct ) and the angle of similarity (0). 


where 


(6.2) 


' a, \ 

/ y 

°oct 

2y/j2 

cos 6 

< 02 1 — ^ 

°oct 

4 7 — 4 

v/3 

cos ($ — 27r/3) 

. CT3 J 

&oct 

K 7 

cos (0 + 27r/3) ^ 


and 


&oct — ^ 


cos30 = -^2 ; °< e < 


(6.3) 


(6.4) 


J 2 can be expressed in terms of a ^ and 6 using equations 6.1, 6.3 and 6.4. 

1 l/m 

I 1 - 


(6.5) 


L b + c ( 373 cos ^ J 

Now, substituting equation 6.5 into equation 6.2 we can express principal stresses in terms 
of just God and 6. 

X / 

COS0 


0 1 


°oct j 

S cr 2 

► = < 

% 

&oct ( 

, CT3 > 


Ooct J 


2 

75 


1 - a(3a oct ) 2mi 
b + c( cos 3 0) 2mi / 3 J 


l/2mi 


cos (0 — 27r/3) 
cos (0 + 2 tt/3) 


(6.6) 


The limits for a oc t can be determined by considering pure hydrostatic tension (7i = 3cr, J 2 = 
J 3 = 0) and compression (7j = -3<r, J 2 = J 3 = 0). Substituting in equation 6.1 gives the 
range as 


-1 

3a 1 / 2mi 


< a < 




(6.7) 
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Now the whole threshold surface is defined by equation 6.1 with a oc t in the range given by 
equation 6.7 and 0 in the range [0,7t/3]. The convexity of F needs to be verified in this 
range. 


The Hessian of F in three dimensional space in matrix form is given by equation 4.32. The 
first and second partial derivatives of F with respect to the invariants can be expressed as 


dF 

dli 

O F 

dJ 2 

dF 

dh 

&F 

dh 2 

d 2 F 

dJ? 

d 2 F 

dtf 


= 2a 1 1 
= b 

— — T- 1/3 
3 3 

= 2 a 

= 0 

- 

9 3 


The derivatives of the stress invariants with respect to the principal stresses can be eas- 
ily found because the invariants are explicitly expressed in terms of the principal stresses 
(equations 3.4). 


The second partial derivatives of F required for the Hessian can now be computed using 

&F_ _ d?F_ dh <Fh dF_ cH _ i <PF chh dJ^ 

doidaj dh 2 do, dcrj dh doidoj dJ 2 2 doi doj 

dF d' 2 J 2 &F_ (hhdJz dF_ d 2 .h 

dJ 2 doidoj dJ' 2 doi doj dJ% doidoj 

where i and j vaxy from 1 to 3. 


The unit normal vector at any point ( 01 , 02 , 03 ) on the threshold surface is 

1 


n = 


v/i& 2 + ii& 2 + © 2 [aai 


dF , dF , dF t 
1 + 7> — 3 a ^ 

0(7 2 OO 3 


(6.9) 


As described in Chapter 4, we need to find two basis vectors in the tangent plane by 
orthonormalizing the unit vectors B , where B = ( Sij-niTij ). U t] in equation 4.31 is the 
matrix containing the two basis vectors in the tangent plane. For convexity, the curvature 
tensor (equation 4.31) must be positive definite at all points on the threshold surface, 
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Figure 6.13: Eigenvalues of the curvature tensor (I 1 J 2 J 3 model) 

which means that the eigenvalues of must be positive. A MATLAB program was written 
to compute these eigenvalues (Appendix C). 

Using the above formulation, a convexity check was performed for all the parameter sets in 
Table 6.8. Very small increments of angle 9 and mean stress <Joct were used in the program 
to exhaustively probe all the threshold surfaces in three dimensional principal stress space. 
Resulting eigenvalues corresponding to each point on the threshold surface were positive 
thereby confirming the positive definiteness of the curvature tensor. Threshold functions 
for parameters corresponding to mi=1.0 in Table 6.8 are convex. The minimum eigenvalues 
of the curvature tensor (for I 1 J 2 J 3 model) at discrete points (using intervals of a 0 ct — 4 MPa 
and 6 = 2 °) are plotted in Figure 6.13. 

The convexity of the function in equation 6.1 was first confirmed for the values of a, b and 
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c (corresponding to m\ =0.6) in Table 6.8. The surface defined by equation 6.1 is observed 
to be identical to that represented by equation 4.23 and hence convexity of the latter is 
established. 


6.3 Biaxial experiments for validation 

The material parameters associated with the flow and evolution («, n, /i, m , (5 , R, H) are 
given in Table 6.7. At this point, we have six sets of threshold surface parameters (a, b , 
c, mi) given in Table 6.8. The focus of this section is the biaxial tests that were used to 
distinguish between these models and thereby choose the proper set of threshold surface 
parameters. 

Generally, in viscoplastic model development, all material parameters are determined using 
uniaxial characterization experiments. Then the model is validated by tests which are very 
different from those used in characterization. These are validation tests that are complicated 
enough to rigorously test the model. 

At this point, we have not conclusively determined the threshold surface parameters. Hence, 
the intent of the biaxial tests is really two fold; (1) to choose the best set of threshold surface 
parameters from Table 6.8 and (2) to validate the model. 

Biaxial tests were chosen for the following reasons: 

1. These tests axe very different from the tensile, compressive and shear tests used in 
characterization. 

2. For non-proportional load paths, the inelastic strain rate vector is not always aligned 
with the loading direction. Hence, such load paths are good tests for the models and 
are expected to differentiate between their predictive capability. 

6.3.1 Test Matrix 

Six biaxial tests (IN4, IN5, IN6, IN2, IN11, IN28), schematically shown in Figure 6.14, were 
performed using strain-controlled axial-torsional loading of aged Inconel 718 tubes at 650° C 
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IN4 

(0.0,0.024) 

100 

(0.02,0.024) 

60 

shear- tensile 

IN5 

(0.018,0.03) 

94 

- 

- 

proportional 

IN6 

(0.03,0.0) 

150 

(0.03,0.024) 

60 

axial-shear 

IN2 

(-0.02,0.0) 

100 

(-0.02,- 

60 

compressive- 




0.024) 


shear 

IN11 

(0.0,0.0042) 

50 

(0.02,0.0125) 

100 

shear- 






proportional 

IN28 

(0.0,-0.0109) 

130 

(-0.015,- 

75 

shear- 




0.0109) 


compressive 
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(experiments performed by C.J. Lissenden at NASA Glenn Research Center). These tests 
were chosen such that they have significantly different load paths in the axial-shear stress 
plane. The loading rates were kept close to the rates used in the characterization tests. The 
strain levels and corresponding times are shown in Table 6.9. The experimental data 
points are shown by circles and squares in Figure 6.15 through 6.28. Though many data 
points were collected only a few are shown to be able to clearly see the comparison with 
predictions of models in Table 6.8. Both stress-time and stress-strain curves are presented. 

6.3.2 Comparison with model predictions 

First, models with mi=1.0 are used to predict stress histories for IN4, IN5, IN6, and IN2 
(Figures 6.15 to 6.24). This section presents and discusses the comparison between model 
predictions and experimental results. 

Shear-tensile loading (IN4) 

Figures 6.15 and 6.16 show the comparison between experiments and model predictions. The 
shear stress predictions during shear loading are identical for I1J2, J2 J3 ■ and I\ J2J3 models. 
This is expected because 7i and J 3 remain zero for shear loading. For the subsequent 
axial loading there is a significant difference in axial stresses predicted by these models. 
Experimental data for axial stress is best predicted by the J2J3 and I1J2J3 models. The J 2 
model over- predicts the axial stresses and so does the I\ J 2 model. The drop in shear stress 
during axial loading is slightly overestimated by all of the models. 

Proportional loading (IN5) 

Axial and shear stresses for the proportional loading are shown in Figures 6.17 and 6.18. 
The J 2 model over-predicts both the axial and shear stresses. Axial stresses are under- 
predicted by I1J2, J2J3, and / 1 J 2 J 3 models and there is no significant difference between 
the predictions. Shear stresses show a moderate difference in model predictions. The I 1 J 2 
model predicts even higher shear stresses than the J 2 model. The I\ J2J3 model predicts 
marginally better than J2 model. The J2J3 model gives the best prediction of shear stresses. 
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Figure 6.17: Proportional loading (IN5): Stress-time response (mi = 1.0) 
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Figure 6.18: Proportional loading (IN5): Stress-strain response (mi = 1.0) 
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Tensile-shear loading (IN6) 


Axial and shear stresses for tensile-shear loading are shown in Figures 6.19 and 6.20. All 
three models {I1J2, J2^3, and /1J2J3) predict almost the same axial and shear stresses. 
Comparison with the experimental data is reasonable but the differentiation between the 
models in terms of their predictive capability is not achieved. Thus, it is difficult to say 
which model gives the best prediction for this experiment. 

There are two possible explanations for the model predictions being almost identical. 

1. Consider variations in the stress invariants /j, J2 and J3 predicted by the J2 model as 
shown in figure 6.21. Axial loading causes all of the invariants Ji, J2 and J3 to increase 
and attain a constant value. The subsequent shear loading causes only a marginal 
increase in J2, while both I\ and J3 drop rapidly. This makes the coefficients a and c 
of the threshold function, which scale I\ and J3 respectively, ineffective. Thus, none 
of the model predictions are significantly different from that of the J2 model. 

2. During shear loading, J3 drops rapidly for all the models and so does I\ for I\ J2 and 
Ji J2J3 models. Ii and J3 become zero approximately at the same time (about 180 
sec) giving a pure shear state (figure 6.21). Thereafter, the material remains in pure 
shear ( J 2 ^ 0, I\ = J3 =0) and the difference between model predictions vanish. 


Compressive-shear loading (IN2) 

Results of the compressive-shear loading are shown in Figures 6.22 and 6.23. Again, all 
models predict nearly the same axial and shear stresses for the same reasons as above for 
IN6. 

During the compressive loading segment of the experiment, the shear stress did not remain 
zero as it should have. A closer look at the data reveals that the shear strains are negligibly 
small (< 10 microstrain). Hence, it is difficult to find a reason for the drift in shear stress. 
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Effective Stress Invariants (MPa) 







This error results in unsatisfactory comparison between predictions and experiment. In the 
event of no experimental error, the model predictions might come close to experiments. 
However, the fact still remains that a good differentiation between the predictive capability 
of the models is not achieved. 

Predictions from the models with rn i=0.6 are shown in Appendix D. There it is observed 
that changing the exponent mi to 0.6 does not change the predictions appreciably. Also, 
assigning different values to mi (other than 1.0 and 0.6) did not change the predictions 
significantly. Thus, further predictions are restricted to models with mi = 1.0. 

6.3.3 Differentiating between models 

Comparing the four biaxial experiments (IN4,IN5,IN6,IN2) it is clear that shear-axial load- 
ing (IN4) is most successful in differentiating between the models and points to a better 
predictive capability of the J2 <^3 and /1 72^3 models. While proportional loading (INS) mod- 
erately differentiates between the models, the load cases starting with axial load (IN6 and 
IN2) are unable to differentiate between the models. 

Figure 6.24 tracks the variations of 7i, J2, and J3 for shear-axial loading (IN4). A clear 
difference between this and Figure 6.21 is that the invariants 7i and J3 increase substantially 
in the second segment (axial) of loading. This makes the coefficients a and c more effective 
and hence we see a large difference in the model predictions. 

A simple way to achieve a substantial increase in I\ and J3 is to start with a shear loading 
and then apply other loading. Based on this idea, two more biaxial experiments (IN11 and 
IN28) were conducted. The load paths for these additional experiments are also shown in 
Figure 6.14. The corresponding test matrix is given in Ihble 6.9. 

Shear-proportional loading (IN11) 

Shear stresses predicted by all models are identical and are in good agreement with the 
experimental data as shown in Figures 6.25 and 6.26 . There is a clear difference between 
the model predictions for axial stress. The J2J3 and I 1 J 2 J 3 model predictions follow the 
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Figure 6.24: Model Predictions: Shear-tensile loading (IN4) 
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wtional loading (IN11): Stress-strain response 
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experimental data closely for most of proportional loading. The I1J2 model over-predicts 
the stress throughout the proportional loading. This further reassurance in the predictive 
capability of the J2J3 and the I\ J2J3 models. 


Shear-compressive loading (IN28) 

Stress predictions for shear-compressive loading are shown in Figures 6.27 and 6 . 28 . Shear 
stresses predicted by all models are not substantially different. The difference in predictions 
is clearly seen for axial stresses. The J2 model over-predicts the compressive stresses and the 
1 1 J 2 model predicts an even higher compressive stress. Predictions using J2 J3 and I\ J2J3 
models show good agreement with the compressive experimental data. 

There is a discrepancy between the experimental values of the shear stress and the model 
predictions during the second loading segment (AB). The shear stress did not drop to almost 
zero as it did in the other experiments (Figure 6 . 15 ). There were no errors noted during the 
experiment and the strain history (shear and axial) is correct. This makes it difficult to say 
why the shear stress did not drop as much as it was expected to. Also, there was a slight 
drift in axial stress in the first loading segment for reasons similar to IN 2 . However, our 
primary interest is in the axial stress during the second loading segment and we disregard 
these discrepancies. This experiment points to the fact that the predictive capability of 
the J2J3 model is better than that of the 7 i J2 model. In this case, the I\ J2J3 model also 
shows a good correlation with experimental data. 


6.4 Summary of model predictions 

Introduction of hydrostatic pressure dependence (through 7 i) in the baseline ( J2) model did 
not improve the model predictions. In the shear-compressive loading (IN 28 ), the 7 i J2 model 
does not predict the axial stress well, even compared to the baseline model. Hence, the J1J2 
model can be ruled out for Inconel 718 . The J2J3 and I\ J2J3 models have consistently 
shown good predictive capability. 

Experiments are in progress at Case Western Reserve University to study the effect of 
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Table 6.10: Material parameters for Incone] 


718 at 650° C 


Parameter 

Value 

k (MPa) 

284.7 

n 

2.0 

H (MPa-s) 

1.7681 xlO 5 

III 

7.0 

P 

3.652 

R (1/s) 

l.Ox 10 -12 

H (MPa) 

1.966 x 10 4 

a (1/MPa 2 ) 

0.0 

b (1/MPa 2 ) 

1.8154 xlO -5 

c (1/MPa 2 ) 

4.5 xlO -6 

mi 

1.0 


hydrostatic pressure on the flow behavior of aged Inconel 718 at room temperature. Initial 
experiments show that the yield and flow in the material is not affected by pressure. These 
experiments suggest that inelastic deformation in Inconel 718 is relatively independent of 
I\ and hence we can choose the J 2 J 3 model over the I 1 J 2 J 3 model. 

All the material parameters for the model are now determined and their values are given 
in Table 6.10. 

6.5 Comments on the method 

Some aspects of the formulation and experimental procedure need to be put in perspective. 

6.5.1 Load Paths 

In three dimensional principal stress space, load paths can be found along which only one of 
the effective invariants changes while the other two remain constant (equation 5.1). Owing 
to the complexity involved with three dimensional testing, two dimensions load paths in 
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Figure 6.29: Principal stress paths with 1\ changing - J 2 and J 3 constant 


axial-shear plane were found along which J 2 is constant and T\ changes (keeping J 3 constant) 
or J 3 changes (keeping I\ constant). Some illustrations of such load paths in effective axial- 
shear stress space were given. However, in the external stress space the corresponding load 
paths are more complicated. This is because the internal stresses evolve during inelastic 
deformation and do not necessarily follow the external stresses. Two load paths in principal 
stress space are shown in Figure 6.29 and Figure 6.30, which correspond to paths CD in 
Figure 5.5 and path AB Figure 5.7 respectively. It may be possible to follow these load paths 
but these tests involve three dimensional loading, which requires complex experimental set 
up. 

In order to avoid the above complexity, an alternative approach was adopted which involves 
simpler experiments. Simple shear experiments helped to determine all of the parameters 
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-400 -1200 

Figure 6.30: Principal stress paths with J3 changing - /1 and J2 constant 


associated with the baseline GVIPS model. The difference is in the determination of the 
additional contribution of I\ (through a), J3 (through c) and exponent mi. Instead of 
trying to determine the contribution of these invariants directly through complicated load 
paths, an indirect procedure involving simple uniaxial and biaxial tests was adopted. A few 
combinations of threshold surface parameters (a, 6, c, mi) were obtained that fit the axial 
test data (tensile and compressive) well. The biaxial tests helped to differentiate between 
the models with respect to their predictive capability. Hence, biaxial experiments were 
useful in identifying the right combination of invariants and also in validating the model. 

The design of experiments in Chapter 5 was helpful in two ways; (1) it established an exper- 
imental procedure to effectively isolate the effect of each invariant on inelastic deformation 
and (2) it provided a methodology for differentiating between the models using biaxial 
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experiments. 


6.5.2 Threshold function 

The choice of the threshold function was the key in formulation. The main advantages of 
this form are addressed below. 

Current form 

There are several advantages of using the generalized polynomial form (equations 4.22 and 
4.23) for the threshold function. 

L The function in equation 4.22 is easily reducible to well known models like Mises and 
Drucker-Prager. 

2. Both functions reduce to the baseline model (equation 4.15) for pure shear loading 
which allows the use of COMPARE for estimation of GVIPS parameters. 

3. The effect of each stress invariant on inelastic deformation can be determined sep- 
arately due to the chosen form of the threshold function. Also, the resulting flow 
and evolutionary laws are mathematically separable into terms associated with each 
invariant. This is a definite numerical advantage in that it allows easy inclusion or ex- 
clusion of a given invariant in the formulation by simply choosing a zero or a nonzero 
value for the corresponding coefficient. 

4. The effectiveness of using the form in equation 4.23 (instead of equation 4.22) was not 
seen in this investigation. This is because the value of the threshold function is not 
significantly different for the two values of m\ (1.0 and 0.6) and the corresponding 
parameters (a, 6, c) in Table 6.8. 

Let us consider an example of axial loading and a I\ model. The values of the 
threshold function are plotted as a function of the ratio a) b (significance of Ti relative 
to J 2 ) for rai=1.0 and mi=0.6 (Figure 6.31). For low values of the ratio (= 0.01 for 
Inconel 718) there is no difference between the threshold value for mi = 1.0 and that 
for mi =0.6. However, as the ratio increases, there is a significant difference between 
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Figure 6.31: Variation of threshold function for axial loading 

the two. Thus, for materials with high a/b ratio (granular or frictional materials), 
changing mi is more effective. For this reason, we retain the form in equation 4.23 
and not revert to equation 4.22 in the formulation. 


Alternative functions 

Some alternative functional forms such as exponential or implicit forms can be considered 
for the threshold function. But these forms will result in flow and evolution equations where 
the invariants are not easily separable. This makes the objective of studying the effect that 
each invariant has on the SD effect more difficult. 

Another possibility is to use a polynomial form with different exponent (equation 4.22 with 
mi t L - 1) for better data fitting capability. Such a function will not reduce to the form in the 
baseline model (equation 4.15) for pure shear loading. Hence the baseline model would have 


NAS A/CR— 2001-210715 


139 




to be altered and the whole optimization software reworked to determine the parameters 
associated with the baseline model. 

6.5.3 Flow law 

Hydrostatic pressure-dependence is represented by the use of I\ in the threshold function. 
An associated flow law for pressure-dependent materials often over-predicts the volume 
expansion during deformation (e.g. Spitzig et al [1975]). However, based on the initial 
experiments on aged Inconel 718 at the Case Western Reserve University, the effect of 
hydrostatic pressure on inelastic deformation appears to be relatively low, resulting in a 
7i- independent material. Hence, the volume expansion is not expected to be over-predicted 
by the associated flow law. Also, the predictions using an associated flow law compare well 
with experiments. Thus, the need for use of a non-associated flow law does not arise for 
Inconel 718. 

6.5.4 Mechanism for SD 

Various deformation mechanisms in Inconel 718 were discussed in Chapter 2 . It appears 
that the particle-dislocation interaction is a more probable cause for the SD than is volume 
expansion. Volume expansion is caused by the retarding effect of hydrostatic pressure on 
the mobile dislocations, with an associated increase in dislocation density. But, since the 
inelastic deformation in Inconel 718 appears to be relatively independent of hydrostatic 
pressure, the volume expansion hypothesis can be ruled out as a cause of SD. Thus, the 
main cause of SD appears to be the interaction between coherency strain fields around 7 " 
precipitates and dislocations. 
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Chapter 7 


Conclusions and Future Work 


This chapter starts by summarizing the steps that were involved in the development of the 
unified viscoplastic model to account for the SD effect. Conclusions that can be drawn 
from this work are presented next. Finally, future work that could be undertaken for the 
betterment of the model is suggested. 


7.1 Summary 

A potential based viscoplasticity model (GVIPS), based on was the starting point of 
the model development. In order to account for the SD effect, the threshold function was 
generalized by including the other two stress invariants (I\ and J3) in its definition. A 
polynomial form was used for the threshold function. Flow and evolution laws to describe 
the material behavior during inelastic deformation were derived from the proposed threshold 
function. 

Experiments were performed on aged Inconel 718 at 650° C to quantify the material param- 
eters in the model. Shear loading was chosen to estimate all the parameters associated with 
the GVIPS model because, for this loading I\ = J3 = 0 and the proposed model reduces to 
the GVIPS model. The parameters associated with I\ and J3 do not affect the determina- 
tion of the other parameters for this loading. Shear tests involved monotonic loading (strain 
and stress control), relaxation, and creep tests with different loading rates and load levels. 


N AS A/CR— 200 1-210715 


141 



GVIPS parameters were optimized to fit the experimental data in shear. The experimental 
data sets were chosen judiciously and the weighting of test data was changed iteratively to 
arrive at the final GVIPS parameters. The goal was to obtain an axial stress-strain response 
intermediate between tension and compression test data. 

Tension and compression tests were conducted up to strains of 2%. By introducing positive 
values for the parameters a (which scales /]) or c (which scales J 3 ) or both, three com- 
binations of invariants were developed. Values of the coefficients were converged upon by 
comparing model predictions with the tension and compression test data. The resulting 
threshold functions had 7i J 2 , J 2 J 3 and I 1 J 2 J 3 forms. Convexity of these threshold func- 
tions in the three dimensional principal stress space was confirmed by checking the resulting 
curvature tensor for positive definiteness. 

An experimental procedure for studying the effect of each stress invariant on the inelastic 
deformation was established in three dimensional principal stress space. This requires fol- 
lowing complicated paths in three dimensional stress space. Simple biaxial tests were chosen 
instead of three dimensional paths for experimental simplicity. Six axial-torsional experi- 
ments were conducted on tubular specimens of aged Inconel 718 at 650°C. Comparison of 
the test data with model predictions indicated that the use of J 3 in the threshold function 
significantly improved the predictive capability of the model especially for non-proportional 
load paths. However, using /j in the formulation did not improve model predictions. In some 
load cases, introduction of I\ widened the gap between experiments and model prediction. 
Also, preliminary results from tension and compression tests conducted under hydrostatic 
pressure indicate that inelastic deformation is relatively independent of pressure, indicating 
I\ -independence. It was therefore established that aged Inconel 718 is a J 2 J 3 material. 


7.2 Conclusions 

The following conclusions can be drawn from the investigation of inelastic behavior of aged 
Inconel 718. 

• A unified viscoplastic model was developed that effectively predicts a strength differ- 
ential effect. 
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• The ability of the model to represent the strength differential effect was demonstrated 
for aged Inconel 718 at 650° C. 

• A general form of the threshold function in terms of the stress invariants I u J 2 , and 

J 3 was proposed and specialized into four trial functions: t/ 2 , I\J 2 t &nd 

The J 2 J 3 threshold function appears to best represent the inelastic response of Inconel 
718, indicating no pressure dependence of flow. 

• The model with a J 2 J 3 threshold function exhibited excellent predictive capability 
under non-proportional axial- torsional loading applied at elevated temperature. 

• A general method to determine the effect that each stress invariant has on inelas- 
tic deformation was established. A current limitation of this method is that three 
dimensional stress paths must be applied. 

• A procedure to determine the convexity of any function in three dimensional stress 
space was implemented to verify convexity of the threshold function. It was established 
that the model developed for aged Inconel 718 satisfies Drucker’s stability criterion. 

• The physical mechanism responsible for the strength differential in Inconel 718 appears 
to be the interaction between coherency strain fields around 7" precipitates and the 
dislocations. 

• The viscoplastic model is very general in that it can be applied to many materials, 
including the ones that do not exhibit a strength differential. The suggested experi- 
mental procedure for characterizing and validating the model is simple and hence can 
be applied to a wide range of materials. 

7.3 Future work 


The viscoplastic model developed in this investigation opens the door for research in related 
areas. 


• More variety of tests can be conducted for model characterization. A test matrix can 
be derived which is best for optimization of the material parameters. 
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• Direct effect of individual stress invariants on the inelastic deformation can be obtained 
by following three dimensional stress paths developed using the proposed threshold 
function. To this end, complicated two dimensional stress paths in the axial-torsional 
stress space with constant J 2 rate, that seem more feasible, can also be attempted. 

• The threshold surface parameters that fit the axial test data well did not fit the 
initial threshold surface data. One way to effectively account for this is to make the 
coefficients of I\ and J 3 (a and c) a function of inelastic deformation. This way, the 
parameters can have initial values that fit the initial threshold surface data, and then 
can evolve such that a good correlation is also obtained with the uniaxial test data. 

• Experimental data for different strain rates can be used in characterizing the model 
so that its predictive capability is improved for a range of loading rates. Similarly, 
if the capability of the baseline model to account for variation in temperature can 
be exploited, a non- isothermal model can be built using this procedure. This would 
additionally involve experiments at different temperature levels. 

• It will be interesting to evaluate the performance of the model using a different class of 
material, e.g. granular materials. The effectiveness of the model would be established 
if it is able to predict the hydrostatic pressure-dependence that is prevalent in these 
materials. 

• Alternative ways of microstructural evaluation of samples deformed under tension and 
compression can be explored to determine the exact mechanism causing SD. 
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Appendix A 


Formulation with Generalized 
Threshold Function 


A.l Flow Law and Evolution Law 


Dissipation Potential 
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Effective Stress Invariants 
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Internal Stress Invariants 
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Derivatives of generalized functions with respect to invariants 
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Derivatives of internal invariants with respect to stress 


l iq a jq 


Derivativ es of generalized functions with respect to stress 
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Substituting from equation A.l 
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Substituting from equation A.l and A. 2 
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Each term in is differentiated by parts and the final expression for the internal compli- 
ance tensor is 
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Simplification for deviatoric models 

For models not involving /1 (hydrostatic pressure-independent material) the 6 k i is inconse- 
quential in the internal constitutive rate equation because 


$kl6tkl = oikk = 0 
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Thus the terms with 5m can be dropped and the internal compliance tensor simplifies to 


Qijkl 
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Stress transformation to principal plane 

Transformation of stresses from the axial-shear plane to the principal plane facilitated the 
numerical inversion of the internal compliance tensor {Qijkl) f° r loadings that start with 
shear. 

A. 2 Contracted Notation 

Various tensors in the formulation have been changed to contracted notation for ease of 
programming. 

6ij= lllOOo] 
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Appendix B 


FORTRAN Program for 
Computing Material Response 


implicit double precision(a-h,o-z) 
common / one / qO ,ql ,q2 ,q3 ,q4 

common / six / xk,xmu,xn,ra, beta, xa,xb,xc,xm,xh, xml, m 
common / eight / ctemp 

dimension sb(6) ,eb(6) ,a(6) ,epl(6) ,dep(6) ,adot (6) , c (6,6) ,acap(6) 

dimension th(6) ,sig(6) ,dev(6) ,dp(6) 

qO = O.dO 

ql = l.dO 

q2 = 2 .dO 

q3 = 3 .dO 

q4 = 4.d0 

pi = q4*atan(ql) 

open (12,file= , nn5.inp , ,status= , old > ) 
open (13,f ile => gnew.dat ’ ,status= > un]mown , ) 

c read current temperature and Poisson ratio 
re ad (12,*) ct emp , xnu 
call GVPROP (XK ,XMU) 
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c read matrix temperature-independent GVIPS parameters 
read(12 ,*) xn,xh,beta,xm,ra 
c read threshold surface parameters 
read (12 , *) xa,xb,xc ,xml 
c material properties (elastic) 

call ELPROPS(emod,alfa,c,xnu) 
gm = emod/(ql + xnu)/q2 
c read no. of load segments and probes 
read(12,*) nstep,ipro 

do 1000, m=l,nstep 

read (12 , *) ttim,ninc , sigi , sigf , taui ,tauf 
timinc = ttim/real (nine) 

tau = taui 
siga = sigi 
deltau = tauf-taui 
delsig = sigf-sigi 
dtau=deltau/real(ninc) 
dsig=delsig/real (nine) 

do 500, k=l,ninc 
c Strain control loading 
eb(l) = siga 
siga = siga + dsig 
eb(2) = tau 
tau = tau + dtau 
ax = emod*(eb(l)-taxe) 
sh = q2*gm* (eb(2)-tshe) 
c transformation to principal stress 
sb(l)=ax/q2 + sqrt (ax**q2/q4 + sh**q2) 
sb(2)=ax/q2 - sqrt (ax**q2/q4 + sh**q2) 
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call GVIPS ( sb , a , adot , dep , c , x j 1 , x j 2 , x j 3 , sig) 
c calculate plastic strains 
do 60, i=l,6 

a(i) = a(i) + adot (i)*timinc 
dp ( i) =dep ( i ) * t iminc 

60 epl(i) = epl(i) + dep(i)*timinc 

axe=q0 

she=q0 

theta=abs(ql/q2*atan(q2*sh/ax) ) 

c transform principal inelastic strain increments to axial and shear components 
axe=ql/q2* (dp ( 1) +dp (2) +abs (dp ( 1 ) -dp (2) ) ♦cos (q2^theta) ) 
she=ql/q2+ (dp ( 1) -dp (2) ) ♦sin(q2^theta) 
c compute total inelastic strain 
taxe=taxe+axe 
tshe=tshe+she 

if (k/1000 . 0-int (k/1000 . 0) . eq. qO) then 
wr i t e ( 1 3 , 999 5 ) k , ct ime ,ax,sh,xjl,xj2,xj3, t axe , t she 
endif 

ctime = ctime + timinc 
500 continue 
1000 continue 
9995 format(i8,9dl2.4) 
stop 
end 

c####################################################################### 
subroutine gvprop(kapa,xmu) 
implicit double precision (a-h,o-z) 
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double precision kapa,kapaO 
common / eight / ctemp 
kapa = 234.7 
xmu = 1.7681d5 
return 
end 

c######################################################################### 

subroutine elprops(emod,alf a,c ,xnu) 
c elastic stiffness tensor 

implicit double precision(a-h,o-z) 
common / one / q0,ql,q2,q3,q4 
common / eight / ctemp 
dimension c(6,6) 
emod = 165360.0 

xlam = xnu*emod/ ( (ql + xnu)*(ql - q2*xnu)) 
xmu = emod/(q2*(ql + xnu)) 


c(l,l) * xlam + q2*xmu 

c(2 ,2) = c(l,l) 

c(3 ,3) = c(l,l) 

c(l ,2) = xlam 

c(l ,3) = c(l ,2) 

c (2 ,3) = c (1 , 2) 

c(2,l) = c(l,2) 

c(3, 1) = c (1 ,3) 

c(3,2) = c(2 ,3) 

c(4,4) = xmu 

c (5 , 5) = xmu 

c(6,6) = xmu 

return 

end 

C############################################################################# 
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subroutine gvips ( st , at , adot , dep , c , x j 1 , x j 2 , x j 3 , sig) 
c Formulation with generalized threshold function 
implicit double precision(a-h,o-z) 
common / one / q0,ql,q2,q3,q4 
common / two / ctemp 

common / six / xk, xmu,xn,ra, beta, xa, xb, xc ,xm,xh, xml ,m 
dimension st (6) ,dep(6) ,dev(6) ,sig(6) ,xl(6,6) , at 0(6) 
dimension acap(6) ,at (6) , adot (6) ,c(6,6),xq(6,6) ,xql(6,6) 
dimension xu(6) ,su(6) ,delta(6) ,xv(6,6) ,xw(6,6) ,yi(6,6) 
dimension pi (6) ,th(6) ,xul (6) ,yil (6 ,6) , thl (6) ,pil (6) ,depi (6) 
dimension indx(6) ,th0(6) ,th2(6) ,pl (6 ,6) ,p2 (6,6) 


if (xa.ne .qO)then 

atO(l) = (Q2*at(l) - at (2) - at(3))/Q3 
at0(2) = (Q2*at(2) - at(l) - at(3))/Q3 
at0(3) - (Q2*at (3) - at (2) - at(l))/Q3 
else 

atO(l)=at (1) 
at0(2)=at (2) 
at0(3)=at (3) 
endif 

at0(4) = at(4) 
at0(5) = at(5) 
at0(6) = at(6) 

xil=at (l)+at (2)+at (3) 
ast=st(l)+st (2) 


C deviatoric stress components 

dev(l) = (q2*st(l) - st(2) 
dev(2) = (q2*st(2) - st(l) 
dev(3) = (q2*st (3) - st(2) 


- st(3))/q3 

- st(3))/q3 

- st (1) )/q3 
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dev (4) = st (4) 
dev(5) = st(5) 
dev(6) = st(6) 

C effective deviatoric stress components 
do 10, i=l,6 

10 sig(i) = dev(i) - atO(i) 

c first effective invariant 

xj l=st ( l)+st (2)+st (3)-xil 

delta(l)=ql 

delta(2)=ql 

delta(3)=ql 

delta(4)=q0 

delta(5)=q0 

delta(6)=q0 

su(l)= sig(l)**q2 + sig(6)**q2 + sig(5)**q2 
su(2)= sig(2)**q2 + sig(4)**q2 + sig(6)**q2 
su(3)= sig(3)**q2 + sig(4)**q2 + sig(5)**q2 
su(4)= sig(6)*sig(5) + sig(2)*sig(4) + sig(3)*sig(4) 
su(5)= sig(l)*sig(5) + sig(4)*sig(6) + sig(3)*sig(5) 
su(6)= sig(l)*sig(6) + sig(2)*sig(6) + sig(4)*sig(5) 

xul(l)= at0(l)**q2 + at0(6)**q2 + at0(5)**q2 
xul(2)= at0(2)**q2 + at0(4)**q2 + at0(6)**q2 
xul(3)= at0(3)**q2 + at0(4)**q2 + at0(5)**q2 
xul(4)= at0(6) *at0(5) + at0(2) *at0(4) + at0(3)*at0(4) 
xul(5)= at0(l)*at0(5) + at0(4)*at0(6) + at0(3)*at0(5) 
xul(6)= at0(l)*at0(6) + at0(2)*at0(6) + at0(4)*at0(5) 

do 12, i=l,6 
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do 12, j=l,6 
xv(i, j)=qO 
xw ( i , j ) =q0 
if (i.eq. j)then 
yi(i, j)=ql 
else 

yi(i, j)=qO 
endif 

12 continue 

c 2nd & 3rd effective invariant 
XI2 = QO 
XJ2 = QO 

xi3 = qo 
XJ3 = QO 

DO 20, 1=1,3 

XI3 = XI3 + Ql* (atO(I) *xul (I)+q2*at0(i+3)*xul (1+3) )/Q3 
20 XJ3 = XJ3 + Ql*(sig(I) *su(I)+q2*sig(I+3)*su(i+3) )/Q3 


DO 21, 1=1,3 

XI2 = XI2 + Ql*(atO(I)*atO(I)/Q2 + at0(I+3)*at0(I+3) ) 
21 XJ2 = XJ2 + Q1*(SIG(I)*SIG(I)/Q2 + SIG(I+3)*SIG(I+3) ) 


sxjl=ql 

sxj3=ql 

if(xjl.lt. qO) sx j l=-ql 
if (x j 3 . It . qO) sx j 3=-ql 

c setting minimum values for internal invariants (Go for stability) 
if (xa.ne . qO . and. xc .ne .q0 . and.ast .ne . q0)then 
cong 1=30.0 
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cong2=900.0 

cong3=27000.0 


else if (xa . ne . qO . and . xc . eq . qO . and . ast . ne . qO) then 
congl=30.0 
cong2=900.0 

elseif (xa . eq . qO . and . xc . ne . qO . and . ast . ne . qO) then 
cong2=900.0 
cong3=27000.0 

elseif (xa . eq . qO . and . xc . eq . qO . or . ast . eq . qO) then 
cong2=900.0 
endif 

if (x j 2 . eq . qO) return 

xek2 = xb*xj2**(xml-ql) 

if (xc .ne .qO .and. ast .ne.qO) then 

xek3 = q2/q3*xc*sxj3*abs(xj3)**(q2/q3*xml-ql) 
else 

xek3 * qO 
endif 

if (xa . ne . qO . and . xc . ne . qO . and . ast .ne . qO) then 
exe=xa*sxj l*abs(xj l)**(2*xml)+xb*xj2**xml+ 

+ xc*sxj3*abs(xj3)**(2*xml/3) 

xekl = q2*xa*sxjl*abs(xjl)**(q2*xml-ql) 

+ -4 . 0/9 . 0*xc*x j 2*sx j 3*abs (xj 3) ** (q2/q3*xml-ql) 

elseif (xc .ne.qO. and. xa.eq.qO .and. ast .ne .qO) then 
exe-xb*xj2**xml+xc*sxj3*abs (xj3)** (2*xml/3) 
xekl =-4 . 0/9 . 0*xc*x j 2*sx j 3*abs (xj 3) ** (q2/q3*xml-ql) 
elseif (xa . ne . qO . and . xc . eq . qO . and . ast . ne . qO) then 
exe=xa*sxj l*abs(xj 1) **(2*xml)+xb*xj2**xml 
xekl = q2*xa*sxjl*abs(xjl)**(q2*xml-ql) 
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elseif (xa. eq.qO. and. xc .eq.qO .or .ast .eq.qO) then 
exe=xb*x j 2**xml 
xekl=qO 
endif 

do 123, i-1,6 

123 pi (i) = xekl*delta(i) + xek2*sig(i) + xek3*su(i) 
c function F 

fcap = (abs(exe))**(ql/xml) - ql 
if (f cap . It . qO) fcap * qO 
flit = fcap**xn 

c inelastic strain rate tensor 
do 30, i=l,6 

30 dep(i) = xk*xk/q2/xmu*f lit* (f cap+ql ) ** (ql-xml) *pi (i) 


if (xa . ne . qO . and . xc . ne . qO . and . ast . ne . qO) then 
if (abs(xil) . It . abs (congl) ) xil=congl 
if (xi2 . It . cong2) xi2=cong2 
if (abs (xi3) . It . abs (cong3) ) xi3=cong3 
elseif (xa . ne . qO . and . xc . eq. qO . and . ast . ne . qO) then 
if (abs(xil) . It .abs (congl ) )xil=congl 
if (xi2.lt .cong2)xi2=cong2 

else if (xa . eq . qO . and . xc . ne . qO . and . ast . ne . qO) then 
if (abs (xi3) . It . abs (cong3) )xi3=cong3 
if (xi2 . It . cong2) xi2=cong2 

else if (xa . eq . qO . and . xc . eq . qO . or . ast . eq . qO) then 
if (xi2 . It . cong2) xi2=cong2 
endif 

sxil=ql 
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sxi3-ql 


if (xil . It .qO)sxil=-ql 
if (xi3 . It . qO) sxi3=-ql 

if (xi2 . gt . qO) then 
if (xc . eq . qO . and . xa . ne . qO) then 
gcap=xa*sxil*abs(xil)**(q2*xml)+xb*xi2**xml 
elseif (xa . eq . qO . and . xc . ne . qO . and . ast . ne . qO) then 
gcap=xb*xi2**xml+xc*sxi3*abs (xi3) ** (q2*xml/q3) 
elseif (xc . ne . qO . and . xa . ne . qO . and . ast . ne . qO) then 
gcap=xa*sxil*abs(xil)** (q2*xml)+xb*xi2**xml 
+ +xc*sxi3*abs (xi3) ** (q2*xml/q3) 

elseif (xc . eq . qO . and . xa . eq. qO . or . ast . eq . qO) then 
gcap=xb*xi2**xml 
endif 
endif 

c function G 

gcap= (abs (gcap) ) ** (ql/xml) 


xik2 = xb*xi2**(xml-ql) 

if ( xc . ne . qO . and . ast . ne . qO ) then 
xik3 = q2/q3*xc*sxi3*abs(xi3)**(q2/q3*xml-ql) 
else 

xik3 = qO 

endif 

if (xc . ne . qO . and . xa . ne . qO . and . ast . ne . qO) then 

xikl = q2*xa*sxil*abs(xil)**(q2*xml-ql) 

+ ~4 . 0/9 . 0*xc*xi2*sxi3*abs (xi3) ** (q2/q3*xml-ql) 

elseif (xc . ne . qO . and . xa . eq . qO . and . ast .ne . qO) then 
xikl = -4.0/9.0*xc*xi2*sxi3*abs(xi3)**(q2/q3*xml-ql) 
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elseif (xa.ne . qO .and. xc . eq.qO . and . ast .ne . qO)then 
xikl = q2*xa*sxil*abs(xil)**(q2*xml-ql) 
elseif (xc . eq.qO. and. xa. eq.qO .or .ast .eq.qO) then 
xikl=qO 

endif 

do 122, i=l,6 

thl (i)=xik2*at0(i)+xik3*xul (i) 

122 th(i) = xikl*delta(i) + xik2*at0(i) + xik3*xul(i) 


xv(l , l)=atO(l) 
xv(l ,6)=at0(6) 
xv(2,2)=at0(2) 
xv(2 ,4)=at0(4) 
xv(3,3)=at0(3) 
xv(3,5)=at0(5) 
xv(4,2)=at0(4) 
xv(4,4)=at0(3) 
xv(5,3)=at0(5) 
xv(5,5)=at0(l) 
xv(6, l)=at0(6) 
xv(6,6)=at0(2) 

xv (1 , l)=atO(l) 
xv(l ,6)=at0(6) 
xw(2 ,2)=at0(2) 
xw(2 ,4)=at0(4) 
xw(3,3)=at0(3) 
xw(3,5)=at0(5) 
xw(4,3)=at0(4) 
xw(4,5)=at0(6) 
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xw(5 , l)=at0(5) 
xw(5,6)=atO(4) 
xw(6,2)=at0(6) 
xw(6 >4)=at0(5) 

if (xc . ne . qO . and . xa . ne * qO . and . ast . ne . qO ) then 
zl = q2*xa*(q2*xml-ql)*sxil*abs(xil)**(q2*xml-q2) 

+ - xb/q3*xi2** (xml-ql) 

+ + 8 . 0*xc/27 . 0* (q2/q3*xml-ql) *sxi3*abs (xi3) ** 

+ (q2/q3*xml-q2) *xi2*xi2 

elseif (xc . ne . qO . and . xa . eq . qO . and . ast . ne . qO) then 
zl=- xb/q3*xi2**(xml-ql) 

+ + 8 . 0*xc/27 . 0* (q2/q3*xml-ql) *sxi3*abs (xi3)** 

+ (q2/q3*xml-q2)*xi2*xi2 

elseif (xa . ne . qO . and . xc * eq. qO . and . ast . ne . qO) then 
zl=q2*xa*(q2*xml-ql)*sxil*abs (xil)**(q2*xml-q2) 
+ -xb/q3*xi2** (xml-ql ) 

endif 


z2 = xb* (xml-ql) *xi2**(xml-q2) 
z3 = xb*xi2** (xml-ql) 

if (xc . ne . qO . and . ast . ne . qO) then 
cnl * sxi3*abs(xi3)**(q2*xml/q3-ql) 
cn2 = sxi3*abs(xi3)**(q2*xml/q3-q2) 
else 
cnl=qO 
cn2=q0 
endif 

cn3 = (ql-xml+beta)/gcap**xml 
c internal compliance tensor 
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if (xa. eq. qO . and . xc . eq.qO . or . ast . eq.qO) then 
do 31, i=l,6 
do 31, j = l,6 

31 xq(i , j )=gcap**beta/xh* 

+ (beta/gcap/xk/xk*atO (i) *atO ( j )+yi (i , j ) ) 

elseif (xc . ne . qO . and. xa. eq.qO . and. ast .ne .qO)then 
do 32, i=l,6 

do 32, j=l,6 

32 xq(i,j) » xk*xk/xh*gcap**(beta+ql-xml)* 

+ (z2*atO(i)*atO( j) + z3*yi(i,j) 

+ + cnl*q2*xc/q3*(xw(i , j)+xv(i,j) 

+ -q2/q3*delta(i) *atO( j ) ) 

+ + cn2*q2*xc/q3* (q2*xml/q3-ql) * (xul (i) *xul ( j ) 

+ -q2/q3*xi2*delta(i) *xul ( j ) ) 

+ + cn3*th(i)*thl(j)) 

elseif (xa.ne .qO . and . xc . eq. qO . and. ast .ne .qO) then 
do 33, i=l,6 

do 33, j =1 , 6 

33 xq(i,j) = xk*xk/xh*gcap**(beta+ql-xml)* 

+ (zl*delta(i)*delta( j ) 

+ + z2*at0(i) *atO( j ) + z3*yi(i,j) 

+ + cn3*th(i) *th( j ) ) 

elseif (xa . ne . qO . and . xc . ne . qO ♦ and . ast . ne . qO) then 
do 34, i=l,6 

do 34, j=l,6 

34 xq(i,j) = xk*xk/xh*gcap**(beta+ql-xml)* 

+ (zl*delta(i)*delta( j ) 

+ + z2*at0(i)*at0(j) + z3*yi(i,j) 

+ + cnl*q2*xc/q3* (xw (i , j ) +xv ( i , j ) -q2/q3*at0 (i) *delta ( j ) 

+ -q2/q3*delt a(i) *atO ( j ) ) 

+ + cn2*q2*xc/q3* (q2*xml/q3-ql ) * (xul ( i) *xul ( j ) - 

+ q2*xi2/q3*xul (i) *delta( j ) -q2/q3*xi2*delta(i) *xul ( j ) ) 
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+ + cn3*th(i) *th( j ) ) 

endif 

c internal stiffness tensor 

call invert (xq, 6,6 ,xl) 

c internal strain rate 

do 70, i=l,6 

70 acap (i ) = (dep ( i) -ra*xk*xk/ xh* 

+ gcap** (ql+beta+xm-xml)*th(i)) 

c convert tensorial shear components to contracted (vector) components 
do 80, i=4,6 

80 acap(i) = q2*acap(i) 

c internal constitutive rate equation 
do 90, i=l,6 
adot(i) = qO 
do 90, j*l,6 

adot(i) = adot(i) + xl(i, j)*acap(j) 

90 continue 

C CONVERT CONTRACTED SHEAR COMPONENTS BACK TO TENSORIAL COMPONENTS 
DO 100, 1=4,6 
100 ACAP (I) = ACAP(I)/Q2 

RETURN 

END 

C ***** * * * * **** ************* 4c *** * * * ** * ** * ** * * ***** * ******** * 4c *********** 4c 
subroutine invert (a, n,m,d) 

c inverts n x n matrix a. result returned as array d 
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c 


m is the dimensioned size of arrays a and d 


implicit double precision(a-h,o-z) 
dimension a(m,m) ,d(m,m) ,w(138,276) 
do 5 i=l,n 
do 5 j = l,n 

w(i, j)=a(i, j) 
if(i.eq.j) then 
w(i,j+n)=l. 
else 

w(i, j+n)=0.0d0 
end if 
5 continue 
do 50 k=l ,n 

do 10 l=l,2*n-k+l 
j=2*n+l-l 

w(k, j)=w(k, j)/w(k,k) 

10 continue 

do 30 i=l,n 

if(i.eq.k) goto 30 
do 20 1=1 ,2*n-k+l 
j=2*n+l-l 

w(i, j)=w(i, j)-v(i,k)*w(k, j) 
20 continue 

30 continue 

50 continue 
do 60 i=l,n 
do 60 j=l,n 

d(i , j )=w(i , j+n) 

60 continue 
return 
end 
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Appendix C 


MATLAB Program for 
Determining Convexity 


m=l .0; 
ij=l; 

for sav=654: -4: -654 
lc=l ; 

for c=2 . 2e-6 
a=2 . 6e-7 ; 
b=1.8154e-5; 
i=l; 

for lam=l : 1 : 27 
ik=i*ij ; 
k=0 ; 

th=pi*lam/180; 

r j2=(l-a* (9*sav~2) ~m) /(b+c/3~m* (2*cos (3*th) ) ~ (2*m/3) ) ; 
rj2=(abs(rj2))~(l/ (2*m)) ; 

7. principal stresses 
sl=sav+2/sqrt (3) *r j2*cos(th) ; 
s2=sav+2/sqrt (3) *r j2*cos(th-2*pi/3) ; 
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s3=sav+2/sqrt (3) *r j2*cos(th+2*pi/3) ; 


7, invariants 
il=sl+s2+s3; 

j 2=1/3* ( (sl-s2) ~2+(s2-s3) ~2+(s3-sl) ~2) ; 

j 3=1/27* (2*sl-s2-s3) * (2*s2-sl-s3)*(2*s3-sl-s2) ; 

snj 3=1 ; 

if j3 < 0 

snj 3=-l ; 

end 

7, derivatives for hessian matrix 
dl=b*j2~(-l + m)*m; 
d3=b*j2~(-2 + m)*(-l + m)*m; 
d5=0; 

if il == 0 
d0=0; 
d6=0; 
else 

d0=a*2*m*il~(2*m-l) ; 

d6=2*a*il~ (-2 + 2*m)*m*(-l + 2*m) ; 

end 

if j3 == 0 
d2=0; 
d4=0; 
else 

d2=(2*c*snj3*(abs(j3))~(-l + (2*m)/3. )*m)/3. ; 
d4=(2*c*snj3*(abs(j3))~(-2 + (2*m)/3.)*(~l + (2*m)/3.)*m)/3. 
end 
djO=l ; 

dj 1=1/3* (2*sl-s2-s3) ; 
dj2=l/3*(2*s2-sl-s3) ; 
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dj3=l/9*(2*sl"2-s2~2-s3~2-2*sl*s2-2*sl*s3+4*s2*s3) ; 
dj4=l/9*(2*s2"2-sl"2-s3"2-2*sl*s2-2*s2*s3+4*sl*s3) ; 
dj 5=2/3; 
dj 6=2/3 ; 

dj 7=2/9* (2*sl-s2-s3) ; 
dj8=2/9*(2*s2-sl-s3) ; 
dj9=-l/3; 

dj 10=2/9* (2*s3-sl-s2) ; 
dj 11=1 ; 

dj 12=1/3* (2*s3-sl-s2) ; 
dj 13=2/3 ; 
djl4=-l/3; 
djl5=-l/3; 

dj 16=1/9* (2*s3"2-sl~2-s2~2-2*sl*s3-2*s2*s3+4*sl*s2) ; 
dj 17=2/9* (2*s3-sl-s2) ; 
dj 18=2/9* (2*s2-sl-s3) ; 
djl9=2/9*(2*sl-s2-s3) ; 

mll=d6 + dl*dj5 + d3*djl*djl + d2*dj7 + d4*dj3*dj3; 
ml2=d6 + dl*dj9 + d3*dj2*djl + d2*djl0+ d4*dj4*dj3; 
m22=d6 + dl*dj6 + d3*dj2*dj2 + d2*dj8 + d4*dj4*dj4; 
ml3=d6 + d3*djl*djl2 + dl*djl4 + d4*dj3*djl6 + d2*djl8; 
m23=d6 + d3*dj2*djl2 + dl*djl5 + d4*dj4*djl6 + d2*djl9; 
m33=d6 + dl*djl3 + d3*djl2"2 + d2*djl7 + d4*djl6"2; 

V, curvature tensor 

x=[mll ml2 ml3 

ml2 m22 m23 

ml3 m23 m33] ; 

ul=dO + dl*djl + d2*dj3; 

u2=d0 + dl*dj2 + d2*dj4; 

u3=d0 + dl*dj 12+ d2*dj 16; 

du=(ul~2+u2"2+u3"2) "0.5; 
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7, unit vector in tangential direction 
u=[ul/du u2/du u3/du ] > ; 
bn^eyeO^u+u’ ; 

7, orthonormalization 
[pm,s,v]=svd(bm) ; 
q=pm ( : ,1:2) ; 

7# curvature in tangential direction 
kur(i)=min(eig(q , *x*q)) ; 
kurr=kur (i) ; 

1 (i)=lam; 

res(ik, l)=sav; 

res(ik,2)=lam; 

res(ik,3)=kurr ; 

if kur(i) < 0.0 

lamc(ik)=lam; 

str(ik)=sav; 

kc(ik)=kurr ; 

cr(lc,ik)=c; 

lc=lc+l ; 

end 

k=k+l ; 

end 

i=i+l ; 

end 

end 
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Stress (MPa) 


Appendix D 


Correlation for Biaxial Loading 



Time (sec) 


Figure D.l: Shear- tensile loading (IN4) (mi = 0.6) 
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(MPa) 



Time (sec) 


Figure D.2: Proportional loading (INS) (mi = 0.6) 
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